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Preface 


The  purpose  of  this  study  was  to  research  and  apply 
current  technology  electronic  data  acquisition  and  tracking 
techniques  to  conventional  munitions  live-fire  testing.  The 
choice  to  apply  Joint  Probability  Data  Association  (JPDA)  in 
a  Kalman  filter  estimator/smoother  permitted  tracking 
multiple  fragments  simultaneously. 

Extensive  Monte  Carlo  simulation  modeling  and  subsequent 
analysis  demonstrated  good  tracking  performance  of  position 
and  velocity  for  non-curving  fragments  and  good  rejection  of 
clutter.  However,  the  simplified  Kalman  filter  first-order 
dynamics  model  did  indicate  difficulty  in  tracking  curving 
fragments,  and  the  overall  results  from  all  simulations 
indicated  poor  acceleration  estimation  performance. 

In  performing  this  research*,  I  am  greatly  indebted  to  my 
entire  thesis  committee;  Capt  Rob  Williams,  thesis  advisor; 
Maj  William  H.  Worsley,  committee  member;  and  Prof  Peter  S. 
Maybeck,  committee  member,  for  their  helpful  encouragement 
and  concerned  feedback  when  it  was  needed  the  most.  A  thank 
you  also  goes  to  Dan  Zambon  of  the  Information  Sciences 
Laboratory  for  tolerating  my  simulation  runs*  bogging  the  VAX 
computers.  My  last  thank  you  goes  to  my  for 

"putting  up  with  it  all"  over  this  past  year. 

Ronald  J.  Beyers 
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The  purpose  of  this  study  was  to  research  and  apply 
current  technology  electronic  data  acquisition  and  tracking 
techniques  to  conventional  munitions  live-fire  testing. 
Previously  applied  high-speed  film  cameras,  celotex  bundles, 
and  associated  technologies  for  munitions  testing  have  proven 
themselves  expensive  in  materials,  labor,  and  time.  Such 
previous  test  methods  cost  upwards  to  $250,000  per  test  blast 
and  require  from  days  to  weeks  to  manually  compile  and  reduce 
collected  blast  data  to  an  analytical  format.  The  specific 
scope  of  this  study  was  to  research  methods  to  electronically 
acquire  and  track  the  position,  velocity,  and  acceleration  of 
multiple  warhead  fragments  as  they  dispersed  from  the  test- 
blast  center.  A  design  specification  for  a  maximum 
trackable  fragment  speed  was  set  at  10,000  ft/sec.  The 
theoretical  application  of  xenon  strobe  illuminated  (2.0 
microsecond  flash  duration) ,  orthogonally  oriented  Charge- 
Coupled  Device  (CCD)  cameras  (in  sets  of  two)  provides  three 
dimensional  image  measurements  at  a  2.0  microsecond  exposure, 
5000  frame/sec  rate.  The  acquired  and  assumed  noisy  fragment 
position  measurements  (recorded  digitally)  are  post-mission 
processed  through  an  Extended  Kalman  Filter  (EKF)  based  Joint 
Probability  Data  Association  (JPDA)  multiple  target/tracker 
state  estimator  followed  by  a  backward  time  Rauch-Tung- 
Striebel  (RTS)  fixed- interval  optimal  smoother.  Strong 
emphasis  was  placed  on  Monte-Carlo  computer  simulation 


testing  of  this  EKF/JPDA/RTS  tracker-smoother  algorithm. 
Representative  trajectories  of  straight,  curving,  and 
crossing  spherical  fragments  at  3000,  6000  and  10000  ft/sec 
were  modeled  and  tracked  with  promising  accuracies  in 
position  and  velocity.  Information  gained  by  these  fragment 
state  variables  allows  the  estimation  of  each  fragment's 
kinetic  energy  and  therefore  lethality.  The  presented 
fragment  data  acquisition  system  was  deemed  realizable  and 
practical  with  existing  technologies,  although  the  CCD  camera 
5000  frame/sec  requirement  was  found  difficult  to  obtain 
reliably.  The  initial  proposed  system  hardware  cost  will  be 
high;  however,  critical  system  components  (such  as  cameras) 
survive  the  test  blast  and  are  continuously  reusable  to  keep 
overall  long-term  costs  down.  In  addition,  the  entire  data 
reduction  process  is  reduced  from  days  or  weeks  to  several 
hours  (overnight)  on  an  autonomous  EKF/JPDA/RTS  computer 
program.  *  V  -  .  '  - 
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JOINT  PROBABILITY  DATA  ASSOCIATION  ON  TRACKING 
MULTIPLE  MUNITIONS  FRAGMENTS 


I.  INTRODUCTION 


I . 1  Problem  Overview 

This  document  presents  research  results  for  an  all- 
electronic  data  acquisition  and  processing  system  for  the 
tracking  of  multiple  airborne  munitions  fragments  for  arena 
munitions  test-blast  applications.  The  problem  is  primarily 
approached  as  a  stochastic  process,  and  employs  recently 
developed  Joint  Probability  Data  Association  (JPDA)  concepts. 
This  research  is  the  first  effort  of  such  a  system  design  and 
therefore  has  a  primarily  theoretical  foundation  combined 
with  numerous  problem  simplifications  to  keep  the  design 
reasonably  tractable. 

I . 1 . a  Problem  Statement.  The  thorough  test  and 
evaluation  of  conventional  munitions  (bombs,  mines,  air-to- 
air  missile  warheads,  etc.)  requires  experimental  detonations 
in  a  test  environment  to  collect  pertinent  functional  and 
lethality  performance  data.  The  3246  Test  Wing  (TW) ,  Eglin 
AFB,  FL  performs  such  test  detonations  as  "arena  tests"  where 
the  test  munition  is  partially  or  fully  surrounded  by  both 
electronic  and  non-electronic  data  collection  devices.  Upon 
detonation,  the  dissemination  of  pertinent  performance  data 
occurs  in  the  micro  and  millisecond  time-frame  and  imposes  a 
very  hostile  environment  toward  the  survivability  of  nearby 


instrumentation  devices.  Past  and  present  arena  test  data 
collection  techniques  have  primarily  employed  passive,  pre¬ 
positioned,  soft  target  arrays  to  catch  a  representative 
portion  of  the  munition  fragment  cloud,  combined  with  minimal 
photo-optical  and  electronic  instrumentation  devices.  Such  a 
data  collection  scheme  is  both  expensive  (in  expended 
materials  and  manhours)  and  of  limited  data  collection 
versatility  and  accuracy  (21) .  This  research  presents  an 
alternative,  all-electronic,  video/ optical  instrumentation 
technology  to  supplement  or  replace  present  data  collection 
systems.  The  intent  is  to  increase  data  quantity  and  quality 
at  lower  total  cost  per  test  in  both  expended  materials  and 
manhours . 


I. l.b  Basic  Approach .  The  general  approach  to  the 
overall  system  design  consists  of  two  parts: 


1)  finding  a  suitable  data  collection/sensing 
technology  for  the  arena  test  mission,  and 

2)  deriving  a  suitable  data  processing  package  that 
applies  stochastic  estimation  techniques  to  extract 
pertinent  fragment  data  from  the  raw  (assumed  noisy) 
data  measurements. 


I.l.b.l  Data  Collection.  The  developed  data 
acquisition  system  employs  multiple  pairs  of  orthogonally 
oriented  Charge  Coupled  Device  (CCD)  video  cameras  to  collect 
three  dimensional  (3D)  position  measurements  of  the  radially 
dispersing  fragments  following  detonation  (see  Appendix  A  for 
geometry) .  For  mathematical  convenience,  this  study  models 
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the  CCD  camera  geometry  such  that,  each  camera  has  a  10.0  x 
10.0  ft  object  plane  at  a  focal  distance  of  20  x  21//2  « 

28.28  ft  in  front  of  the  camera.  Appendix  A  gives  a  full 
description  of  the  camera  and  arena  geometry  model.  The  use 
of  CCD  cameras  is  derived  from  a  thorough  study  of  sensor 
alternatives  that  is  discussed  in  chapter  II.  Based  on  arena 
test  reports  that  indicate  maximum  fragment  speeds  near  8000 
ft/sec  (16;  18),  and  an  additional  20%  design  margin,  this 
system  is  designed  for  a  maximum  expected  fragment  speed  of 
10,000  ft/sec.  This  maximum  speed  specification  defined  the 
need  for  a  high  measurement  (frame)  rate  with  a 
corresponding  short  exposure  time  for  each  frame  to  "freeze" 
motion.  These  characteristics  of  high  frame  rate  combined 
with  the  short  exposure  time  per  frame  are  deemed  feasible 
through  the  use  of  synchronous  xenon  stroboscopic 
illumination  with  the  test  performed  at  night.  Several  CCD 
camera  pairs  and  associated  strobe-lamp  arrays  are  placed 
throughout  the  test  arena  to  provide  the  necessary  volume 
coverage  necessary  to  collect  sufficient  fragment  data,  as 
determined  by  munitions  analysts.  The  collected  image  data 
is  recorded  (magnetic  tape  or  optical  disk)  during  the  test 
blast  for  post-mission  reduction  and  analysis.  Figure  1.1 
presents  this  data  acquisition  scheme.  In  this  figure,  the 
timing  and  control  computer  controls  the  timing  and 
synchronization  of  camera  frame  rate,  strobe  lamp  triggers, 
and  other  instrumentation  timing.  The  collected  analog 


camera  images  are  digitized  (A/D  converters)  and  routed  to  an 
applicable  high-speed  digital  recorder  (several  magnetic  tape 
drives  or  laser  optical  disk  system) .  All  data  is 
simultaneously  recorded  with  the  master  clock  and  appropriate 
control  status  flags. 


Cam  1A  Cam  IB  Cam  2 A  Cam  2B .  Cam  nA  Cam  nB 


Figure  1.1.  Master  Data  Acquisition  Block  Diagram 
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I . 1 .b. 2  Data  Reduction/Processina .  The  recorded 
image  data  is  played  back  and  processed  through  a  fragment 
state  estimator  package  presented  in  this  document  consisting 
of  an  Extended  Kalman  Filter  (EKF)  operating  together  with  a 
Joint  Probability  Data  Association  (JPDA)  algorithm  (see 
Figure  1.2).  The  derived  forward-time  state  estimates  are 
then  reverse-time  processed  through  a  Rauch-Tung-Striebel 
(RTS)  fixed-interval  smoother  to  optimize  the  fragment  state 
estimates  for  all  time  from  all  information  available  within 
the  recorded  measurement  data  associated  with  the  test-blast 
event.  Estimated  states  include  the  time  varying  three- 
dimensional  (3D)  position,  velocity  and  acceleration  for 
each  fragment  observed  by  the  cameras. 


Figure  1.2.  Post-Mission  Data  Processing  Block  Diagram 


Knowledge  of  the  above  states  allows  further  analysis  to 
determine  each  fragment's  time  varying  kinetic  energy  and 
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associated  lethality.  Such  lethality  determination  is  not 
the  purpose  of  this  research,  only  the  initial  design  and 
simulation  of  the  instrumentation  system. 


1.2  Scope 

This  research  addresses  a  number  of  the  important 
aspects  regarding  the  data  acquisition  and  tracking  of 
multiple,  high-speed  targets  (fragments).  It  is  emphasized 
that,  at  the  time  of  this  research,  no  previous  work  had  been 
found  or  reported  regarding  this  task.  Based  on  information 
extracted  from  previous  arena  munitions  test  reports  and 
correspondence  with  the  3246  Test  Wing,  elementary  level 
design  specifications  are  developed  to  define  the  overall 
technical  problem  (e.g.  the  10,000  ft/sec  spec.).  From  these 
specifications,  a  suitable  acquisition/ sensing  technology  is 
selected,  followed  by  the  development  of  the  previously 
mentioned  EKF/JPDA/RTS  Tracker/ Smoother . 

1.3  Assumptions 

The  following  assumptions  are  made,  primarily  to 
simplify  this  initial  investigation  of  the  proposed  approach 
and  to  illustrate  the  important  basic  concepts  involved. 
Another  reason  for  making  these  assumptions  is  the  lack  of 
data  needed  to  support  assuming  otherwise.  Many  of  the 
assumptions  are  presented  during  the  development  portion  of 
this  document;  however,  some  of  the  main  primary  simplifying 
assumptions  are  stated  here. 


1)  Fragments  generally  follow  a  straight-line 
trajectory  from  the  arena  center  (location  of  test 
munition)  radially  outward  with  minimal  lateral 
accelerations  normal  to  the  fragment's  velocity 
vector. 

2)  The  only  accelerations  acting  on  a  fragment  (post 
detonation)  are  those  due  to  aerodynamic  drag  and 
gravity. 

3)  Fragments  are  modelled  as  spheres  consisting  of 
uniform  material  density. 

4)  Each  fragment  is  traveling  in  the  transonic  to 
hypersonic  (0.9  to  4.0  Mach  number)  speed  range,  and 
their  coefficients  of  drag,  C  ,  are  set  to  1.0 
based  on  explicit  test  data  for  spherical 
projectiles  presented  in  Figure  1.3,  below: 


(22:742) 

Figure  1.3.  Zero-Lift  Drag  Coefficient  for  Various  Bodies; 

where  M  is  the  relative  speed  in  Mach. 


5. )  All  sensors  (cameras,  etc)  are  perfectly  aligned 

and  calibrated  with  no  vibration  induced  from  the 
munition  detonation.  This  significantly  simplifies 
the  measurement  noise  modeling  requirements. 

6. )  Each  frame  of  collected  image  data  from  the  CCD 

cameras  has  been  pre-processed  to  determine: 

a)  Each  observed  fragment's  image  centroid 

coordinates,  (x_,y_). 

c  c 

b)  Each  observed  fragment's  diameter,  d. 

c)  The  total  integer  number  of  fragments  present  in 
each  camera's  image. 

These  three  measurement  variables  are  then  supplied 
to  the  EKF/JPDA/RTS  package  for  processing. 

7. )  All  noises,  unless  otherwise  stated,  are  considered 

white,  Gaussian,  stationary,  and  mutually 
independent . 

From  the  above  set  of  assumptions,  a  simplified  yet 
reasonably  realistic  set  of  models  and  solutions  may  be 
derived. 


1.4  Chapter  Overviews 

1.4 .a  Background .  The  background  covers  a  general 
review  of  present  arena  test  technologies  and  their  relation 
to  the  problem  statement.  A  research  review  for  various 
sensor  technologies  is  given  to  describe  the  CCD  camera 
selection  as  the  applicable  sensor.  This  is  followed  by  a 
brief  review  of  Kalman  filter  basics  and  the  EKF  structure  by 
applicable  equations.  This  leads  into  a  background  and 
structural  description  of  the  JPDA  algorithm,  including  a 
general  example.  The  last  section  gives  a  brief  overview  of 
the  Rauch-Tung-Striebel  fixed-interval  optimal  smoother 
structure. 
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1.4 .b  Development.  This  chapter  explains  the  applied 
system  reference  frames,  followed  by  the  derived  fragment 
dynamics  truth  model  description.  A  detailed  equation  layout 
follows  of  the  applied  Singer  approach  used  in  the  EKF 
structure.  A  breakout  of  the  EKF  nonlinear  measurement  model 
is  then  given,  accompanied  by  a  measurement  ambiguity 
example.  The  last  section  describes  in  detail  the  applied 
EKF  initialization  procedure  developed  in  this  study. 


1.4 . c  Simulation.  This  chapter  gives  a  general 
description  of  the  developed  simulation  programs  with  their 
respective  features.  In  addition,  this  chapter  lists  the 
actual  numerical  values  and  constants  used  in  the  equations 
developed  in  previous  chapters.  A  detailed  description  is 
included  of  the  applied  pseudo-random  noise  corruption 
process.  The  chapter  ends  with  a  representative  example 
output  plot  generated  by  the  simulation  package. 

1.4. d  Analysis.  This  chapter  explains  in  detail  the 
specific  simulation  details,  nuances,  and  findings.  This 
chapter  also  contains  the  analysis  of  special  "extra  runs" 
used  to  investigate  peculiarities  noted  in  the  EKF/ JPDA/RTS 
algorithm. 

1.4. e  Conclusions  and  Recommendations .  Gives  an  overall 
listing  of  good  and  bad  performance  characteristics  gathered 
in  the  research.  This  is  followed  by  a  list  of  interest 
items  that  need  further  attention  in  regards  to  more 
simulation  runs  or  system  design  modifications. 


This  background  chapter  contains  three  main  parts.  The 
first  presents  the  present  arena  test  data  acquisition  system 
in  use  by  the  3246  Test  Wing.  The  second  deals  with  a 
research  and  applicability  study  of  various  sensor 
technologies  for  the  arena  test  application.  The  final  part 
discusses  Kalman  filters;  putting  emphasis  on  past  work 
regarding  the  JPDA  algorithm. 

II. 1  Past  Testing  Methods 

11. 1.  a  Test  Profile.  Isolated  portions  of  the  Eglin  AFB 
test  range  are  reserved  for  the  conduct  of  arena  blast  test 
missions  (14;  15:225).  Munitions  tested  include 
omnidirectional  and  shaped  charge  munitions  that  range  in 
size  from  small  antipersonnel  (AP)  mines  and  shells  to  2000 
pound  general  purpose  (GP)  bombs.  Additionally,  for  a  single 
given  munition  type,  several  tests  may  require  differing  sets 
of  data  collection  under  differing  test  conditions.  As 
discussed  in  Chapter  I,  fragment  velocities  commonly  approach 
8000  ft/sec.  Therefore,  the  overall  conventional  munitions 
arena  test  mission  requires  an  adaptable  data  collection 
scheme  that  will  accommodate  a  large  variation  of  munitions 
type,  test  conditions,  and  range  of  fragment  velocities. 

11.1. b  Present  System.  The  present  primary  data 


acquisition  system  in  use  by  the  3246  Test  Wing  for  its  arena 
testing  is  the  Vulnerability  and  Lethality  Testing  System 


(VALTS) ,  a  computerized  system  that  employs  light  sensors, 
velocity  screens,  and  pressure  sensors  (21) .  Supplementing 
VALTS  is  the  use  of  fire-resistant  mineral  fiberboard  target 
bundles,  high-speed  color  film  motion  picture  cameras,  and  a 
standard  closed  circuit  TV  camera  for  safety  purposes. 

Light  sensors  of  the  VALTS  are  used  to  detect  the  first 
light  that  defines  the  initiation  of  case  breakup  (time  tQ) . 
This  initial  t  becomes  starting  time,  to  which,  all 
subsequent  measurements  and  observations  are  referenced  to. 
Timing  information  is  maintained  by  generating  and  recording 
an  Inter-Range  Instrumentation  Group  (IRIG)  time-code  with 
all  VALTS  measurements. 

A  VALTS  velocity  (Dahlgren)  screen  is  normally  placed 
immediately  in  front  of  each  fiberboard  bundle  to  detect  each 
fragment's  exact  arrival  time  at  the  bundle  (see  Figure  2.1). 
The  Dahlgren  screen  senses  the  fragment  arrival  when  the 
fragment  provides  an  electrically  conductive  path  between  the 
sprayed  aluminum  and  aluminum  foil  layers  as  it  passes 
through  the  structure.  (21:Sec  III).  The  resulting  "pulse” 
is  recorded  with  the  IRIG  time-code  of  the  pulse,  t  ,  on 

tr 

magnetic  tape  for  post-mission  analysis.  The  average 
velocity  for  each  detected  fragment  may  then  be  determined  by 
taking  the  known  distance  from  the  test  munition  to  the 
velocity  screen,  r  ,  and  dividing  it  by  the  elapsed  time 
since  t  (21:9): 

Vavg  =  rs  /  (tp  -  V  W 
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Figure  2.1.  Exploded  View  of  Typical  Dahlgren  Velocity  Screen 

Piezoelectric  pressure  sensors  of  VALTS  are  normally 
placed  in  various  locations  throughout  the  arena  to  detect 
relative  blast  pressures  and  to  record  their  time  of 
occurrence  with  the  IRIG  time-code  (21:24).  The  fiberboard 
target  bundles  may  vary  in  size  for  each  test.  Typically, 
the  bundles  are  eight  feet  high,  four  feet  wide,  and  stacked 
four  feet  thick.  They  are  positioned  to  intercept  a 
representative  segment  of  the  fragmentation  cloud  as  it 
disperses  from  the  detonated  munition.  The  bundles  are  then 
post-mission  disassembled  one  layer  at  a  time  and  the 
coordinates,  weight,  size,  and  composition  of  each  embedded 
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« 

fragment  are  tabulated  for  later  analysis  (16:79-202;  18:21- 
28;  20 : APPX  B) .  • 

Lastly,  each  arena  test  normally  includes  at  least  one 
high-speed  motion  picture  camera  operating  at  up  to  8000 
frames  per  second  optically  recording  the  test  munition  ^ 

detonation  and  its  effect  on  the  surrounding  arena 
environment  (16:4;  17:18;  18:6;  19:16;  20:3).  Example  arena 
set-ups  are  shown  in  Figures  2 . 2  and  2.3:  .1 


(20:5) 

Figure  2.2.  Example  Arena  Set-up,  General  Overhead  View 
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(18:10) 

Figure  2.3.  Detail  of  Arena  Set-up,  with  Example  Dimensions 

II. 2  Review  of  Candidate  Sensor  Technologies 

II. 2. a  Radar/Lidar  Imaging.  Various  radar  technologies, 
including  standard  wavelength  (less  than  35  GHz)  and 
millimeter  wavelength  (40-300  GHz),  did  not  meet  the  required 
spatial  resolution  required  to  individually  track  such  a 
large  number  of  closely  co-located,  high-speed  fragments. 

This  resolution  shortfall  is  primarily  due  to  oversized 
beamwidth  caused  by  limited  hardware  capabilities.  Laser 
radar  can  theoretically  satisfy  the  spatial  resolution 
requirement  with  its  inherently  narrow  beamwidth;  however, 
the  required  number  of  illumination  lasers  and  their 
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associated  pointing  and  scanning  slew- rates,  based  on  the 
10,000  ft/sec  tracking  requirement,  rendered  such  a  system 
unrealistically  complicated  and  expensive  (7;  40:553-566). 
Overall,  any  sort  of  radar  fails  to  meet  the  design  criteria 
in  either  spatial  resolution  or  dynamic  tracking  capability 
(bandwidth) . 

II . 2 .b  Optical  Imaging.  The  implementation  of  optical 
imaging  type  sensing  systems  can  theoretically  satisfy  the 
necessary  three  dimensional  spatial  resolution  and  tracking 
requirements  as  long  as  the  frame  rate  is  fast  enough  to 
prevent  any  appreciable  motion-induced  smearing  of  the 
fragment  image  (3;  25:313-318).  Such  imaging  devices  may  be 
classified  into  two  general  classifications:  Photo- 
electronic  image  devices  (camera  tubes) ,  and  solid-state 
(chip)  imaging  devices. 

II.2.b.l  Electronic  Image  Devices  (Tubes) .  The 
most  common  photo-electronic  image  device  applied  to 
instrumentation  applications  is  the  vidicon  tube,  although 
other  tube  type  imagers  exist  such  as  the  image  isocon,  image 
orthicon,  plumbicon,  and  secondary  emission  conduction  (SEC) . 
All  these  camera  tubes  have  existed  for  decades  and  basically 
operate  similarly,  in  that  an  optical  image  is  focused  on  a 
photosensitive  layered  photocathode  (image  target)  which 
generates  a  corresponding  video  signal  by  scanning  a 
relatively  low  current,  low  velocity  electron  beam  across  the 
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image  target  (3;  26;  27).  The  dominant  drawbacks  of  any 
direct  optical  input  camera  tube  used  for  recording  high¬ 
speed  events  are  beam-discharge  and  photoconductive  lag  of 
the  target's  photosensitive  layer  —  causing  image  smearing 
(36?  41:24-29).  The  more  dominant  beam-discharge  lag  is  due 
to  the  incomplete  neutralization  of  the  entire  image  target 
surface  by  a  single  scan  of  the  electron  beam,  as 
graphically  presented  in  plot  (c) ,  of  Figure  2.4  (36;  41). 
Photoconductive  lag  is  defined  where  the  layer's  charge  does 
not  change  instantaneously  with  a  corresponding  instantaneous 
change  in  light  input,  as  shown  graphically  in  plot  (b) ,  of 
Figure  2.4.  The  important  facts  to  realize  from  these  plots 
in  Figure  2.4  are  that  they  represent  a  0.02  second  frame 
period  (50  frames/sec)  and  that  the  arena  test  application 
requires  a  much  shorter  frame  period  of  200  Msec  (5,000 
frames/sec) .  Therefore,  the  slow,  exponential  decay  rates 
displayed  in  this  figure  would  dominate  the  camera  tube's 
operation  and  render  its  output  a  uselessly  smeared  image  in 
such  a  high  frame-rate  application.  A  commonly  applied 
technique  that  is  actually  a  spill-over  from  high-speed 
photography  is  to  apply  high-speed  shuttering  between  the 
image  and  the  target  to  "gate''  the  incoming  light  (12).  Such 
a  technique  shortens  the  target  exposure  time:  effectively 
"freezing"  the  moving  image,  with  the  only  additional 
requirement  being  increased  illumination  levels  to  compensate 
for  the  shortened  photo-optical  integration  period  of  the 
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T  ■  0.02  sec 
(50  Hz  Frame 
Rate) 


(■)  The  target  potential  7l  before  the  moment  of  stabilization  and  the  target 
stabilisation  potential  Fa  aa  a  function  of  the  signal  output  J  V.  (b)  The  target  potential 
of  a  picture  element  as  a  function  of  time  if  the  light  intensity  is  switched  from  white  to  a 
dark  grey,  at  t  =  0.  (c)  The  signal  output  current  at  successive  scans. 


Figure  2.4. 


(36:242) 

Characteristic  Plots  of  Beam-Discharge  Lag 
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image  target.  Mechanical  rotating  disk  type  shutters  can 
effectively  shorten  the  target  exposure  time  down  to  one 
millisecond  (msec) .  Looking  at  the  worse  case  tracking 
velocity  requirement  of  10 , 000  ft/sec  combined  with  a  10  ft 
wide  field  of  view  (FOV) ,  a  1.0  msec  gate  period  would  result 
in  a  10  ft  streak  image  of  that  fragment  which  is  the  full 
FOV.  Therefore,  the  gate  period  must  be  shortened 
significantly  for  effective  ’’freezing"  of  the  10,000  ft/sec 
fragment.  Shortening  the  gate  period  to  2.0  /isec  results  in 
a  streak  of  0.02  ft  («0.25  inches),  which  is  acceptable. 

This  2.0  nsec  gate  period  can  be  achieved  by  using 
electronically  gated  image  tubes  or  stroboscopic  illumination 
techniques  (3;  12;  33).  However,  electronically  gated  image 
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tubes  are  undesirably  bulky,  expensive,  and  fragile.  They 
also  require  awkward  support  circuits  and  voltages?  an 
undesirable  choice  for  the  arena  test  environment. 
Stroboscopic  lighting  techniques  perfected  for  high-speed 
photography  can  be  directly  applied  to  electronic  imaging. 

The  only  unusual  requirement  would  be  that  the  arena  tests 
would  have  to  be  conducted  at  night  since  there  must  be  no 
ambient  light  between  strobe  flashes  (12:1;  34;  42). 

Once  the  gate  period  problem  is  resolved  and  the 
exposure  time  is  restricted  to  the  desired  2.0  nsec,  the 
induced  image  must  be  read  off  the  image  tube's  target  within 
the  remaining  198  nsec  before  the  next  strobe  flash. 

Recalling  the  previous  discussion  regarding  target  lag,  to 
read  and  neutralize  the  target  and  then  to  be  ready  for  the 
next  flash  within  198  nsec  is  virtually  impossible  even  with 
the  fastest  responding  image  tube  (plumbicon)  (12:5).  An 
analogy  to  this  dilemma  would  be  taking  a  snap-shot  with  a 
standard  film  camera  but  not  being  able  to  remove  the  film 
before  another  exposure  is  taken  on  top  of  it. 

This  leaves  the  final  conclusion  that  tube  type  imaging 
devices  cannot  satisfy  the  required  200  nsec  frame  rate. 

This  is  true  even  though  their  exposure  gate  period  can  be 
reduced  to  the  desired  2.0  nsec  period. 

II . 2 .b. 2  Solid-State  Imagers .  Solid-state  imaging 
devices  include  charge-coupled  devices  (CCD's)  and  charge- 
injection  devices  (CID's)  (4;  10;  30;  38).  Sensitivity  and 
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resolution  characteristics  have  steadily  improved  over  the 
years  (30) .  Both  solid-state  technologies  are  similar  in 
that  the  input  image  is  focused  on  a  silicon  planar  array 
(chip)  of  photosites.  The  photon  bombardment  on  each  of 
these  photosites  induces  an  accumulation  of  minority 
carriers  proportional  to  the  light  intensity  on  each 
photosite.  These  accumulated  charges  are  periodically 
transferred  out  of  each  photosite  and  eventually  off  the  chip 
by  appropriate  logic  switching  signals.  The  off-loaded 
signals  are  then  amplified  and  mixed  with  appropriate 
synchronization  signals  to  produce  a  usable  video  signal. 

Charge  coupled  devices  differ  from  CIDs  in  the  actual 
mechanism  used  to  transfer  the  accumulated  photocharges.  In 
CCDs,  the  individual  photocharges  are  transferred  ("coupled") 
from  one  charge  well  to  another  towards  the  perimeter  of  the 
chip  (similar  to  a  "bucket  brigade"  of  firemen  passing 
buckets  of  water) ,  where  an  output  circuit  then  converts  the 
charges  to  minute  currents  for  off-chip  processing  (30) .  The 
CID  differs  in  that  the  photocharges  are  directly 
transferred  ("injected")  to  the  chip's  substrate  for  direct 
output  (10) .  Both  devices  and  their  associated  charge 
readout  mechanisms  have  their  tradeoffs  in  operational 
flexibility  and  performance,  dependent  on  the  proposed 
imaging  application. 

In  a  study  to  compare  the  relative  spectral 
sensitivities  and  resolution  characteristics  of  the  vidicon 
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tube,  CCD,  and  CID,  the  solid-state  devices  outperformed  the 
vidicon  in  all  categories  except  where  the  vidicon  resolution 
(28.4  lines/mm)  outperformed  the  CCD  resolution  (23.0 
lines/mm) .  The  CID  performed  better  under  low  light 
conditions,  while  the  CCD  performed  better  with  scene 
contrast  values  of  50  per  cent  or  greater  (29:44).  The  arena 
testing  environment  will  unlikely  be  a  low  light 
classification,  since  additional  illumination  may  be  provided 
if  warranted  by  the  system  design.  This  fact  favors  the  CCD 
over  the  CID.  Another  important  drawback  with  the  CID  is  its 
larger  inherent  readout  capacitance  (10:192)  which  limits  the 
device  readout  speed  to  well  below  the  required  198  fjiS 
specification.  This  effectively  eliminates  the  CID  as  a 
candidate  technology. 

Current  applications  of  CCDs  are  widespread,  especially 
in  military  applications  as  surveillance,  guidance,  and 
target  identification  (6;  9;  28;  31).  In  addition,  CCDs  have 
a  strong  potential  commercial  market,  which  has  motivated 
commercial  research  and  development  (38) .  Present  CCD 
resolution  densities  of  1280  x  970  pixels  (photosites)  are 
available  on  the  commercial  market  at  production  cost 
(2:122) . 

Using  the  CCD  instead  of  an  imaging  tube  leaves  the  same 
problems  of  requiring  a  2.0  Msec  optical  gating  period  and 
198  m sec  to  read  out  the  acquired  image  from  the  chip  before 
the  next  image  must  be  taken.  CCDs  can  and  have  been 
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coupled  with  electronically  gated  image  intensifier  tubes  as 
standard  image  tubes  have.  However,  as  previously  pointed 
out,  gated  image  tubes  are  undesirable  for  the  arena  test 
application.  Again,  using  stroboscopic  lighting  as  the 
gating  mechanism  appears  the  most  technically  and 
economically  feasible.  Although  the  majority  of  CCD  imager 
chips  have  been  designed  for  slow  (0.02  sec)  frame  readout, 
there  have  been  chip  architectures  and  techniques  developed 
to  allow  full  frame  readout  well  within  the  198  /xsec  period 
(1) .  Thermocouple  cooling  of  the  CCD  chip  to  approximately 
-40*  C  significantly  improves  the  CCD's  signal-to-noise 
ratio,  thereby  raising  its  optical  sensitivity  (4:142). 

Many  of  the  above  described  characteristics  of  CCDs  make  them 
the  desirable  choice  for  the  arena  test  application. 

11 . 3  Final  Sensor  Selection 

Radar,  lidar,  tube  imaging,  and  solid-state  imaging 
technologies  were  investigated  as  possible  solutions  to  the 
problem  of  tracking  10,000  ft/sec  fragments.  After  analyzing 
each  technology,  the  following  conclusions  were  made: 

1.  Use  thermocouple  cooled  CCD  cameras  as  the  image 
sensing  device  operating  at  5000  frames/sec. 

2.  Use  stroboscopic  illumination  to  facilitate  2.0  nsec 
optical  gating  to  "freeze"  fragment  motion. 

11. 4  Kalman  Filters 

The  reader  is  assumed  to  have  a  thorough  knowledge  of 
Kalman  filtering  concepts.  This  section  is  only  a  review  of 
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the  basic  concepts  and  equations.  For  a  more  in-depth 
review,  see  (33) . 

The  original  Kalman  filter  formulation  can  only  be  used 
for  linear  systems,  while  the  arena  test  application  of 
interest  is  partially  nonlinear.  The  linear  form  has  been 
revised  through  perturbation  techniques  to  handle  nonlinear 
systems.  This  revision  brought  about  two  formats:  1)  the 
linearized  Kalman  filter  and  2)  the  extended  Kalman  filter. 
The  linearized  filter  assumes  a  nominal  trajectory  over  the 
entire  estimation  space  of  interest,  while  the  extended 
filter  revises  the  trajectory  declaration  based  on  the  latest 


state  estimate.  The  linearized  filter  is  simpler  to 
implement,  making  it  attractable  for  on-line  computer 
application.  However,  its  largest  drawback  is  its  tendency 
to  break  track  if  the  measurements  deviate  significantly  from 
the  precalculated  nominal.  The  extended  filter  overcomes 
this  off-nominal  sensitivity  (within  reasonable  limits)  at 
the  cost  of  being  more  computationally  burdensome.  For  the 
application  at  hand,  on-line  computational  capability  is  not 
required.  In  addition,  although  a  series  of  nominal 
trajectories  for  the  arena  test  environment  (i.e.  there  are 


more  than  one)  may  be  pre-generated  and  stored,  the  large 
quantity  of  "typical  nominals"  would  yield  a  linearized 
filter  that  would  quickly  surpass  the  computational  load  of 
the  extended  filter's  implementation.  The  result  is  to  apply 
the  extended  filter  format  including  a  linear  dynamics, 
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single-input  (gravity) ,  and  nonlinear  measurement  model  as 
the  following  equation  outline  illustrates. 


I 


For  a  given  linear  stochastic  single-input  system,  its 
dynamics  model  (in  continuous  form)  may  be  defined  as: 


x(t)  =  F(t) x(t)  +  b(t) u (t)  +  G(t)w(t)  (2.1) 

where 

3c(t)  =  the  n-dimensional  vector  of  time-varying  system 
state  variables, 

x(t)  =  the  time  derivative  of  x(t) , 

F(t)  =  the  n-by-n  dimensional  system  dynamics  or 
"plant"  matrix, 

b(t)  =  the  n-dimensional  control  input  gain  vector, 

u(t)  =  the  scalar  control  input  (gravity), 

£(t)  =  the  n-by-s  dimensional  noise  distribution  and 
gain  matrix, 

w(t)  =  the  s-dimensional  vector  of  mutually 

independent  white  Gaussian  noises  of  expected 
value: 


E{  w(t)  }  =  0  (zero-mean),  (2.2) 

and  of  strength  fl(t) ,  where 

fl(t)  =  the  s-by-s  dimensional  w(t)  strength  matrix 
defined  by: 

E{  w(t)  w(t')  }  =  fl(t)  &(  t  -  t’),  (2.3) 

and 

E{  •  )  =  the  expected  value  of  {  •  ) 

5 ( • )  =  the  delta  function  at  time  ( • ) . 
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The  previous  system,  Eq.(2.1),  rewritten  into  discrete 
propagation  equation  form  becomes: 


2E(tk)  4(tk,tk-l)  -^k-l5  + 
+  Sd(tk-1}  Hd(tk-1^ 


rk 

J~(t' 

L-tk-l 


r)  fe(T)dr 


u<tk-l) 

(2.4) 


where  #  (t^,  is  the  state  transition  matrix  associated 

with  F(t)  :  the  solution  to  which  is  described  by  (32:40): 


In  Eq. (2.4),  5d(*)  may  be  dropped  by  setting  it  equal  to  the 
identity  matrix  without  loss  of  generality: 


fidC)  “  I 

allowing  Eq.(2.4)  to  simplify  to  a  final  state  propagation 
equation  as: 


x(tR)  -  *(t]C'tk-1)  ^k-l*  +  fed(tk-l}  u^tk-l^  +  -d(tk-l) 

(2.5) 


where 


«.(*,•)  =  discrete-time,  zero-mean,  white,  Gaussian 

noise  sequence  independent  of  X(t^) ,  and  has 
the  following  covariance  kernel: 


■{^(V^tj))  = 


0 


where 


for  k  *  j 
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fld(tk)  =  J  *(tk,r)  G(r)  fi(r)  G~(r)  i‘(tk,r)  dr  (2.6 
t,_  , 


The  Kalman  filter  state  estimate  propagation  equation  based 
on  this  model  is: 


where 


*<V  =  *<Wi>  ‘iVii  +  Wi>  u(tk-i)  <2-7 


tk  denotes  the  propagatated  time  just  prior  to  update 


tk-1  denotes  the  post-update  time  of  the  previous 
recursion  cycle. 


Meanwhile,  the  discrete-time  propagation  form  for  the 
filter's  error  covariance  becomes: 


*  *<Wi>  *<<-!>  *T<Wi>  + 


+  j*(tk,r)  O(r)  9(r)  fiT(r)  lT(tk,l)  dr  (2.8 


E(tk)  =  *«£-!>  +  VW  <2-9 


For  the  stochastic  nonlinear  discrete  measurement  model: 


z(tR)  =  h[x(tk),tk]  +  Y(tR) 


(2.10 


where 


lit*#*]  *  an  m-dimension  vector  of  functions  of  2£(tk) 
and  time, 

Y(*#*)  =  an  m-dimensional  vector  of  mutually 

independent,  zero-mean,  white  Gaussian 
discrete  noise  sequences,  independent  of 
x(tk)  and  H(*r*)#  that  has  the  following 
covariance  kernel: 


E{v(tk)  v^tj)} 


B(fck) 

o 


for  k  /  j 


(2.11) 


The  tine-variant  state-space  to  measurement-space 
transformation  matrix,  H[x(tk),tk),  that  is  needed  within 
the  extended  Kalman  filter  is  defined  as: 


H[x(tk),tk]  = 


bfilSSn'V 


“  4<tK> 


(2.12) 


The  update  equations  for  the  residual,  r(tk) ,  appear  as: 


r(tk)  =  zR  -  h[x(tk),tk] 
where 


(2.13) 


z.  is  the  actual  measurement  data  (a  realization 
of  Eq. (2.10) ) 

h[x(t~)  ,tk]  is  the  best  prediction  of  the  measurement  value 
before  it  arrives  at  sample  time  tk  . 


The  residual  covariance,  S(tk)  ;  Kalman  qain,  K(tk)  * 
updated  state  estimate  vector,  x(tk) ;  and  updated  estimate 


covariance  matrix,  P(tk)  are  now  evaluated  as: 

S(tk)  =  H[*(tk),tk)  P(tk)  HT[i(tk),tk]  +  g(tk)  (2.14) 
S(tk)  =  P(t~)  HT[x(tk),tk]  S_1(tk)  (2.15) 
*(tk)  =  *(t“)  +  K(tk)  r(tk)  (2.16) 
£(tk)  =  P(tk)  -  K(tk)  H[k(tk),tk]  £(t")  (2.17) 


where 

and 


superscript  denotes  pre-update  values 
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+  superscript  denotes  post-update  values 
Further  details  of  the  KF  development  including  actual 
numerical  assignment  methods  are  given  in  Chapter  III. 

II. 5  Joint  Probability  Data  Association  f JPDA1  History 

The  JPDA  algorithm  is  actually  an  outgrowth  from  the 
Probabilistic  Data  Association  (PDA)  algorithm  developed  by 
Bar-Shalom  and  Tse  (5) .  The  PDA  and  JPDA  algorithms  are 
intended  for  applications  involving  one  or  a  combination  of: 

1)  Multiple  simultaneous  measurements 

2)  Multiple  targets  (including  false  targets  or 
"clutter”) 

3)  Multiple  sensors 

Such  situations  define  a  track/data/measurement  association 
problem  that  PDA  and  JPDA  are  intended  to  solve. 

There  are  two  ways  to  approach  the  track  association 
problem,  from  a  measurement  oriented  approach  or  a  target 
oriented  approach.  The  measurement  oriented  approach, 
represented  in  Figure  2.5,  associates  each  measurement  with 
all  possible  prior  target  tracks,  new  tracks,  or  clutter. 

From  this  association,  a  hypothesis  tree  is  made  to  determine 
the  most  likely  association  of  the  new  measurements  with 
previously  determined  tracks.  As  the  tracks  are  updated, 
appropriate  weights  are  assigned  to  compute  the  conditional 
probabilities  (37)  . 

The  target  oriented  approach,  represented  in  Figure  2.6, 
uses  the  established  target  tracks  to  gain  statistical 
information  about  the  propagated  ex\  "ted  fragment  position. 
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This  information  is  then  combined  "a  posteriori”  with  the 


Figure  2.5.  Target  Weighted  Tracking  Method 


Figure  2.6.  Measurement  Weighted  Tracking  Method 
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The  previous  articles  explained  various  applications  for 
measurement  oriented  and  target  oriented  approaches.  In 
addition,  the  work  accomplished  by  Purvis  (35),  which  used 
the  target  oriented  approach,  was  studied  for  comparison  to 
the  arena  test  application.  The  arena  test  application,  with 
its  multiple  target  and  potentially  high  clutter  environment 
was  compared  to  these  published  examples.  The  target 
oriented  approach  was  deemed  most  applicable  to  the  fragment 
tracking  application  since  the  literature  indicated  it  to  be 
better  at  rejecting  clutter.  Therefore,  the  target  oriented 
JPDA  algorithm  was  chosen  for  this  design.  The  combination 
of  work  by  Bar-Shalom,  Chang,  Fortmann,  and  Tse  (5;  11;  24) 
set  the  foundation  for  the  EKF/JPDA/RTS  reduction  package 
developed  in  this  research. 

II . 6  Joint  Probability  Data  Association  (JPDA)  Specifics 
The  originally  devised  PDA  scheme  is  primarily  intended  to 
track  a  single  target  in  a  cluttered  environment.  Although 
PDA  has  been  applied  to  tracking  scattered  multiple  targets, 
its  reliability  becomes  questionable  for  target  clustering  or 
for  the  worse  case  of  track  crossing  (5) .  The  original  PDA 
discarded  measurements  as  clutter  if  they  did  not  associate 
with  a  target.  Therefore  the  PDA  scheme  was  modified  to 
include  probability  terms  for  any  measurement  not  associated 
with  a  certain  target,  then  that  measurement  might  be 
associated  with  a  different  target.  The  result  is  a  set  of 
joint  probabilities  across  the  entire  measurement  and  track 
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event  space  including  clutter.  Referring  to  Figure  2.7 ,  the 
resulting  JPDA  probabilities,  p,  are  incorporated  as 
weighting  terms  in  the  update  process.  This  is  accomplished 
by  the  generation  of  modified  residuals,  j£;  residual 
covariances,  S;  and  filter  covariances,  £;  based  on  all 
feasibly  validated  permutations  and  combinations  (events)  of 
measurement  to  target/track  association  including  clutter. 
The  result  is  a  weighted  residual  sum  of  all  residuals  that 
might  correspond  to  each  target.  The  residual  covariance 
becomes  the  key  player  in  measurement-to-target  association. 
The  larger  the  residual  covariance  for  a  particular 
measurement/target  pair  the  lower  its  probability  term  and 
contribution  to  the  weighted  sum  of  residuals.  This  process 
will  become  clear  in  the  following  section  that  details  the 
equation  steps. 


Figure  2.7.  JPDA  System  Block  Diagram 


(35) 
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The  implementation  of  the  JPDA  follows;  a  detailed 
derivation  is  found  in  (24)  .  Key  definitions  include; 


n  =  the  number  of  targets 

M  =  the  number  of  measurements 

m  =  the  dimension  of  the  measurement  vector 

i  =  the  target  index  (i  =  0,  1,  2...n) 

j  =  the  measurement  index  (j  =  1,  2,  3...M) 

X  =  index  of  each  feasible  event;  defined  as  no 

more  than  one  measurement  originates  from  each 
target  and  each  measurement  has  only  one  origin 


Tj(X) 


=  measurement  association  indicator; 


Tj(X) 


1  if  measurement  j  is  associated 
with  any  target  for  event  X 

0  otherwise 


£^(X)  =  target  detection  indicator; 


«i(X) 


1  if  any  measurement  is  associated 
with  target  i  for  event  X 

0  otherwise 


P { X 1 Z  }  =  the  probability  that  event 
all  measurements  Z  where 


X  is  true  given 
Z  =  Z[t0,tk]. 


Note  that  the  previous  target  index  counter,  i,  begins 
at  "0"  to  indicate  the  "no  target"  or  clutter  condition.  The 
JPDA  algorithm  begins  by  generating  the  above  defined  0 
validation  matrix.  At  update  time  t^,  the  Q  is  generated 
using  a  "g-sigma”  ellipsoid  test  (24)  beginning  with  the 
filter  residual  similar  to  Eg. (2.13): 

rj(tk)  =  Zj(tk)  -  h^xft")  ,t‘]  (2.17) 

where  in  this  case 

r^t.  )  =  the  residual  corresponding  to  the  association 
^  of  measurement  j  and  the  predicted 

measurement  for  target  i  at  time  t~. 

and  the  g-sigma  test  now  equates  as: 


I4  j 


(tk)] 


[  sx(tk) 


r1 


< 


(2.18) 


where 


[  s1(tk)  ] 


and 


9 


=  the  inverse  of  the  residual  covariance  matrix 
for  target  i  as  defined  by  eq. (2.13) 

=  the  m-dimensional  validation  gate  volume 
defined  to  have  a  chi-squared  probability 
distribution. 


An  example  pair  of  two-dimension  (m  =  2)  validation 


gates  that  would  be  generated  by  Eq.(2.18)  and  the 
corresponding  Q  matrix  generated  from  the  pictured  case  of 
four  candidate  measurements  (M  =  4),  through  and 


appears  in  Figure 


two  current  targets  (n  =  2) ,  and  h2» 

2.8.  For  each  element  ur.^  of  ft,  the  g-sigma  test, 

Eq.(2.18)  is  performed  to  determine  which  measurement/target 

pairs  are  valid  (set  that  or.^  element  equal  to  one  if  true; 

zero  if  false).  A  "good"  measurement,  j,  will  form  a  small 

residual  when  tested  against  the  correct  target,  i;  producing 

a  small  product  in  Eq.(2.18)  and  "falling  inside"  the  volume 
.  .  .  2 

limit  defined  by  g  .  Also  note  m  the  pictured  Q 
validation  matrix  that  the  first  column,  i  =  0,  corresponds 
to  each  measurement,  z^  not  associated  to  any  target,  h^>Q 
therefore  making  it  clutter.  This  first  column  is  always 
filled  with  ones  since  there  is  always  some  probability  that 
each  measurement  may  be  clutter.  Based  on  the  previously 
stated  definition  of  a  feasible  event, 


Figure  2.8.  Example  Validation  Gates  with  Associated 
ft  Validation  Matrix 


the  Q  matrix  roust  be  broken  down  into  individual  Q(X) 
event  matrices  made  up  of  elements  drj^(X).  This  event 
generation  process  is  governed  (24:176)  by  the  following 
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rules: 


1)  Scan  n  by  rows  and  pick  one  unit  per  row  for  Q(X) 
(i.e.  there  can  be  only  one  origin  for  a 
measurement. 

2)  Only  one  unit  from  each  column  i  >  1  can  be  taken 
(i.e.  at  most  one  measurement  could  have  originated 
from  a  target) .  The  number  of  units  from  column 

i  =  0  is  not  restricted  because  any  measurement  can 
be  clutter. 

Although  at  first  glance  the  previous  rules  for 
generating  the  individual  Q(X)  matrices  may  appear  trivial, 
much  thought  is  required  to  develop  a  permutation/combina¬ 
tion  event  algorithm.  The  devised  algorithm  must  ensure  all 
measurement-to-target  permutations  as  well  as  all  possible 
clutter  combinations  are  obtained  for  any  given  Q  matrix 
without  inadvertent  repeats  of  events.  Appendix  B  gives  the 
mathematical  formulation  for  determining  the  maximum  number 
(worse-case)  of  possible  events  given  M  measurements,  n 
targets,  and  a  fully  populated  (no  zero  elements)  n  matrix. 

For  the  previous  example  Q  matrix  (Figure  2.8),  if  the 
n  matrix  were  fully  populated  with  ones,  then  Eq. (B.7)  of 
Appendix  B,  indicates  there  would  be  21  possible  events. 
However,  due  to  the  validation  gating  process  setting  some 
0  elements  to  zero,  this  number  has  been  reduced  to  12 
possible  Q(X)  event  matrices,  thereby  "trimming"  the  size 
of  the  hypothesis  tree.  The  resulting  twelve  Q(X)  event 
matrices  defined  by  the  above  rules  appear  as: 
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0  0  1 

10  0 

10  0 

0(io)  = 

10  0 

Q(ii)  = 

10  0 

0(12)  = 

10  0 

0  10 

0  10 

10  0 

From  each  event  matrix  n(X),  the  measurement 
association  indicator,  r ^ (X) ,  is  defined  equal  to  one  if  any 
element  (other  than  the  first  column)  in  measurement  row  j 
of  fi(X)  equals  one;  or  it  is  set  to  zero  otherwise. 
Likewise,  the  target  detection  indicator,  6^(X),  is  set  to 
one  if  any  element  of  target  column  i  (except  for  i  ■  0)  o 


Q(X)  equals  one;  or  it  is  set  to  zero  otherwise. 

The  next  step  is  to  calculate  the  event  probabilities, 
P{X|z  },  for  each  event  where: 


« 


I 


I 


exp  -0.5  [f^(tk)]T  tS1(tk)]"1  rj(tk^ 


1 


where 

=  the  probability  of  detection  for  target  i 
'  (assumed  equal  for  all  targets  in  this  study  to 
simplify  application) , 

C  =  the  density  of  false  measurements,  assumed 
Poisson  distributed  by  (24) , 

c  =  normal iz at ' ''n  constant, 

0(X)  =  the  number  of  false  measurements  for  event  X  found 
by: 


(  2  *  )*/2  JS1 (tk) | 1/2 


<  1  -  PD  > 


(2.19) 


i:6^  =  0 


C*(X) 

P( x | zk }  =  - 


c 


i: S .  -  1 


0(X) 


T.  t  i  -  t.(x)  ] 

j=i  J 


(2.20) 


Next  calculate  the  association  probabilities,  py  for  each 
event  from: 


«  Z1  P(X|zk)  ur..(X) 
J  X  J 


2  —  1,  2,  3,...  M 

(2.21) 

i=0,  1,  2,...  n 


while  the  associated  probability  that  target  i  goes 
undetected,  is: 
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I 


i 

0 


1  - 


i  =  0,  1, . . .  n 


(2.22) 


From  the  preceding  terms  the  combined  weighted  residual  is: 


(2.23) 


which  is  then  applied  into  target  i's  account  of  Eq. (2.16) 
as  the  state  estimate  update.  Meanwhile,  for  target  i's 
filter  covariance  update,  Eq.(2.17)  is  modified  to  account 
for  clutter  and  measurement  origin  uncertainty  to  re-equate 
as: 


PX(t£)  =  Px(tk)  -  (1  -  /3j)  Ei(tk)  Sx(tk)  [K1(tk)]T  +  JP' 1  (tk) 

(2.24) 

where  the  superscript  i  indicates  target  i's  account  of 
Eqs . (2.14)  and  (2.15).  For  the  correction  matrix,  E'X(t.): 


P’^V  =  KX(tk) 


-  rx(tk)  [r1(tk)]T 


[K1(tk)]T 


(2.25) 


This  concludes  the  JPDA  structure. 

II. 7  Rauch-Tuna-Striebel  (RTS)  Fixed-Interval  Optimal 
Smoother 

This  section  gives  a  brief  summary  of  the  Rauch-Tung 
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Striebel  (RTS)  fixed-interval  optimal  smoother  structure. 
For  a  detailed  derivation  and  analysis  refer  to  (33) .  The 
RTS  smoother  has  a  continuous-time  and  a  discrete-time 
format.  The  latter  of  the  two  formats  was  the  applicable 
choice  for  the  presented  EKF/JPDA/RTS  tracker/smoother 
system. 

Computationally,  the  RTS  algorithm  entails  performing 
the  forward  EKF  computations  and  storing  the  x(t~)  »  E(t^)  , 
x(t£)  and  P(t*)  values  for  all  time  over  the  interval 

[to^fl* 

The  final  forward-time  updated  estimate,  *(t^)  is 
applied  as  the  starting  boundary  condition  and  the  smoothed 
estimate  is  generated  backward  in  time  as  (33): 


x(tf|tf)  «  x(tj)  (2.26) 

X(tkltf)  -*(t+>  +  4(tk)[*(tk+1|t£)  -*ct£+1f]  (2  2?) 

where  the  smoother  gain  matrix  is  given  by: 

A(tk)  =  £(t+>  lT(tk+1,tk)  E~1(tk+1)  (2.28) 

Note  that  in  Eq.(2.28),  the  time  indices  for  the  #(*,*) 

matrix  seem  to  indicate  a  forward-time  transition  matrix. 

However,  because  of  the  adjoint  nature  for  this  system,  the 
T 

transpose,  #  (•,•)  reverses  the  propagation  direction 
(8:276).  Also  note  special  care  must  be  exercised  in 
software  implementations  that  P(»)  is  checked  to  be  non¬ 
zero  before  the  matrix  inversion  is  attempted. 


38 


For  the  reverse  propagation  covariance  generation,  the 
final  time  updated  forward  filter  covariance,  E(t^) ,  becomes 
the  initial  value  so: 


P(tf|tf)  =  P(t+) 


(2.29) 


and 

E(tkltf)  =  £(<>  +  4(tk)[£(tk+1|t£)  -  P(tk+1)]AT(tk) 


(2.30) 


The  implementation  of  the  RTS  smoother  is  quite  simple 
because  of  its  reverse-time  recursive  structure  as  the 
preceding  equations  indicate. 


II. 8  Summary 

This  chapter  has  detailed  the  arena  test  background  and 
its  current  technologies.  This  was  followed  by  a  research 
review  that  compared  various  sensing  technologies  and 
resulted  in  the  selection  of  the  CCD  camera  for  this 
application.  The  remainder  of  the  chapter  reviewed  the 
mathematical  structures  of  the  extended  Kalman  filter,  joint 
probabilistic  data  association  algorithm,  and  Rauch-Tung- 
Striebel  optimal  smoother. 
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III.  DEVELOPMENT 

The  design  and  validation  of  the  developed  system  is 
determined  by  the  fabrication  of  a  computer  simulation  for 
dynamic  "exercise”  of  the  applied  algorithm.  Before  the 
simulation  can  be  accomplished,  parameterized  models  had  to 
be  derived  that  would  reasonably  represent  typical  arena  test 
profiles.  In  addition,  the  resulting  profiles  must  remain 
within  the  limits  of  the  simplifying  assumptions  in  Chapter 
I.  This  chapter  outlines  the  applicable  reference  frames  and 
details  the  development  of  1}  the  fragment  trajectory  truth 
model,  and  2)  both  the  EKF  dynamics  and  measurement  models. 

The  real-world  applicability  of  the  designed  package 
can  only  be  judged  against  the  fidelity  of  the  truth  model 
that  is  used  as  the  real-world  reference.  However,  no  all- 
inclusive  closed  form  solution  for  hypersonic  projectiles 
passing  through  a  fluid  could  be  found  other  than  the 
standard  zero-lift  drag  equation: 

D ( t)  =  0.5  CD  S  O  |v(t) |2  (3.1) 

where  ' 

D(t)  =  the  force  due  to  aerodynamic  drag  (lbf) 

CD  =  the  coefficient  of  drag  for  the  fragment  (unitless) 

2 

S  =  the  cross-sectional  area  of  the  fragment  (ft  ) 

2 

Note:  S  =  0.25  it  d  where  d  =  the  spherical  fragment 

diameter  (ft) 

.  3 

pa  =  the  atmospheric  density  (lbm/ft  ) 

|y(t) |  =  the  magnitude  of  the  relative  airspeed  velocity 
vector  (ft/sec) 
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t  =  time  (seconds) 

This  equation  is  selected  as  the  applied  truth  model  for 
the  fragments. 

III. 1  Reference  Frames 

A  geometry  model  consisting  of  four  reference  frames, 
associated  coordinate  axes,  and  transformation  equations  are 
devised  for  the  arena  to  facilitate  mathematical  tractability 
in  the  algorithm.  For  the  proposed  arena  set-up,  each 
camera  pair  defines  three  of  its  own  reference  frames,  a  two- 
dimensional  (2D)  image  frame  for  each  camera,  and  a  three- 
dimensional  (3D)  intermediate  frame  formed  at  the 
intersection  of  the  two  image  frames.  The  collective  set  of 
all  3D  intermediate  frames  are  then  transformed  into  the 
world  coordinate  frame  located  at  the  arena  center.  A 
complete  description  of  the  arena  geometric  model  is 
contained  in  Appendix  A. 

This  study  only  models  one  camera  pair  and  therefore 
includes  the  following  four  reference  frames: 

1)  Two  camera  image  plane  frames,  one  for  each  camera  as 
(x^y^  and  (x2,y2). 

2)  The  intermediate  coordinate  frame,  (Xj,  Yj,  Zj) , 
located  at  the  orthogonal  intersection  of  the  two 
camera  image  planes. 

3)  The  world  coordinate  frame,  (Xw,  Yw,  Zw) ,  located  at 
the  center  of  the  arena . 
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The  EKF/JPDA/RTS  package  operates  in  the  intermediate 
reference  frame  to  reduce  computational  loading,  and  the 
resulting  estimates  would  normally  be  transformed  to  the 
world  coordinate  frame.  To  minimize  excessive 
transformations  between  the  world  and  intermediate  frame, 
this  study  performed  all  error  analysis  in  the  intermediate 
frame.  Figures  A. 3,  A. 4,  and  A. 5  in  Appendix  A,  illustrate 
the  above  reference  frames. 


Ill . 2  Fragment  Dynamics 

III . 2 . a  Truth  Model .  Based  on  the  simplifying 


assumptions  outlined  in  the  previous  chapter,  a  free-body 


representation  of  a  typical  fragment,  including  the  forces 
acting  upon  it  as  it  travels  through  the  atmosphere  in  free 


flight , 


appears  as: 


localJ 

VERTICAL 

VELOCITY 

VECTOR 

DRAG  INDUCED 
ACCELERATION  4 

^-ARBITRARY 

ANGLE 

> 

'GRAVITATIONAL 

ACCELERATION 

Figure  3.1.  Free-body  Diagram  of  Typical  Fragment  with 
Arbitrarily  Pointed  Velocity  Vector. 


The  magnitude  of  the  force  applied  to  the  fragment  due  to 
aerodynamic  drag,  Eg. (3.1),  is  repeated  as: 
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D(t)  =  0.5  CD  S  ^>a  |v(t)  |2 

From  this  drag  force,  the  associated  drag  acceleration 
upon  the  fragment  is  found  by  dividing  the  force  D(t)  by 
the  fragment  mass.  For  the  case  of  uniform  density  fragments 
(as  assumed)  the  overall  mass  is  expressed  as  a  function  of 
diameter  as: 

m-  =  jt  d3  /  6  (3.2) 

where  1 

m^.  =  the  fragment  mass  (lbm) 

3 

=  the  fragment  material  density  (lbm/ft  ) 
d  =  the  fragment  diameter  (ft) 


Equations  (3.1)  and  (3.2)  may  now  be  combined  to  form  an 

equation  giving  the  drag-induced  acceleration  magnitude 
2 

aD(t)  (ft/sec  )  as  a  function  of  velocity  magnitude  (speed), 
diameter  and  time. 

'3  Pa  CD  1  *2 

a  (t)  =  - -  (3.3) 

Vf d 

The  negative  sign  indicates  deceleration.  For  notational 
convenience,  we  can  define 


K1  " 


-V a  CD 


Pf 


yielding  Eq.(3.3)  more  simply  as 


(3.3a) 


= 


K2  I  V(t) 


(3.3b) 


From  the  above  nonlinear,  speed  dependent,  fragment 
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acceleration  model,  fragment  trajectories  are  generated 


starting  with  an  initial  world  frame  position  E(tQ)  and 
velocity  v(to) .  This  is  accomplished  by  the  standard  two 
integration  steps  of  acceleration  and  velocity  to  generate 
new  acceleration,  velocity,  and  position  vector  values  as 
functions  of  time: 


Figure  3.2.  Fragment  Trajectory  State  Generation. 

Ill . 3  Filter  Dynamics  Model 

The  dynamic  model  must  basically  follow  a  free-falling 
mass  (fragment)  in  the  intermediate  reference  frame  where  the 
high  initial  velocity  and  its  associated  drag  acceleration 
clearly  dominates  any  gravitational  acceleration  component. 
The  EKF  must  estimate  three  states,  position,  velocity  and 
acceleration  in  each  of  the  Xj,  Yj,  and  Zj  directions  (total 
of  nine  states)  for  the  fragment  of  interest.  The  filter 
structure  begins  as  three  separate,  three-state  filters 
augmented  together  to  form  a  single  nine-state  EKF. 
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Originally,  three  separate,  uncoupled,  three-state  filters 
were  designed,  but  had  to  be  coupled  together  as  one  large 
filter  because  of  the  JPDA  algorithm  requires  the  existence 
on  one  residual  covariance  matrix.  The  Singer  dynamics 
model  (23;  39)  is  selected  because  of  its  good  performance  in 
tracking  radar  algorithms  demonstrated  in  (23)  for  tracking 
mild-maneuvering  targets.  Because  the  fragments  are  assumed 
to  be  following  a  nearly  straight-line  trajectory,  they  may 
be  classified  as  mild-maneuvering  targets.  Such  an 
assumption  fits  well  with  the  first-order  Markov  acceleration 
model  included  in  the  Singer  approach.  The  important  fact  to 
keep  in  mind  during  the  following  outline  of  the  Singer  model 
is  that  all  the  necessary  numbers  will  be  provided  by  the 
initialization  routine  discussed  later  in  this  chapter.  The 
dynamics  equations  following  Eq.(2.1)  with  G(tk)  =  X  are; 
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I  «• 

I 


0 

0 

0 

0 

0 

wx 

0 

0 

0 

•  u(tk)  + 

0 

0 

wy 

0 

0 

1 

0 

0 

w 

Z_ 

where  _ 

u(t.  )  =  -32.1741  ft/sec  (gravity) 

(constant  for  all  t^  ) 

r  ,  t  ,  and  t  =  the  first-order  acceleration  time 
x  y  z  constants,  set  different  for  each 

respective  direction  as  determined  by 
the  initialization  procedure  to  yet  to 
be  discussed 

w  ,  w  and  w  =  mutually  independent,  zero-mean,  white, 
y  Gaussian  noise,  each  with  a  covariance 

kernel  as  defined  in  Eq.(2.6)  and 
detailed  below  in  Eqs.(3.8)  and  (3.9). 


Converting  to  the  discrete-time  format  as  outlined  in 
Eq. (2.6)  results  in  the  following  (nine-by-nine)  forward-time 
state  transition  matrix  (13) : 


• 

*x  * 

0 

•  0 

±(T)  = 

0  . 

*y 

•  0 

0  . 

• 

0 

•  iz_ 

• 

where  *(T)  is  parameterized 

by 

rx'  Ty' 

as: 


(3.5) 
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Ti[T/Ti  "  1  +  exp(-T/ri)3 
rL[  1  -  exp(-T/r i) ] 
exp(-T/ri) 


(3.6) 


where 


i  =  x,  y(  or  z  and 


T  -  -  Vi 


(3.7) 


=  200.0  juSec  (constant  sample  period) 


For  the  discrete  bd(tk)  input  gain  vector,  all  elements  are 

2 

zero  except  for  elements  b?  and  bg  which  are  T  /2  and  T 
respectively.  This  bd(tk)  constant  f°r  all 

The  covariance  matrix  or  noise  strength  for 

wd(tk)  is  block  diagonal  symmetric  (nine-by-nine)  for  the 
Singer  approach  and  defined  as: 


W  = 

•  *  •  « 
0 


(3.8) 


where 


2  a2. 
ai 


T5/20 

T4/8 

T3/6 

T4/8 

T3/3 

T2/2 

T3/6 

2 

T  /  2 

T 

x,  y, 

or  z 

2  2  ,2 
“ax'  °ay  and  °a2 


(3.9) 


the  variances  or  mean-squared  values 
of  the  fragment  accelerations  for 
each  respective  direction  determined 
normally  by  "a  priori"  knowledge  or 


The  initial  filter  estimate  covariance  matrix  P(t  ) 

O' 

for  the  Singer  approach  is  also  block  diagonal  symmetric 
(nine-by-nine)  defined  as: 


E(t0) 


where 


Pl.l  Pl,2 

P2 , 1  P2 , 2  P2 , 3 

0  P3,2  P3 , 3 


and  each  indicated  element  of  Eg. (3.11)  is  expressed  as: 


Plfl  =  amm(V 


Pl,2  =  P2 , 1  "  <WV  /  T 


2  aL(t  )  o2.  (t  )  rj  r 


-  2  ETERM  - 


ETERM 


(V  Ti  T  T 

- - - -  ETERM  +  - 


where 


P3,3  “  °ai<V 


(3.16) 


i  =  x,  y,  or  z 

2 

a  (t  )  =  the  position  measurement  corruption  noise 
°  variance  at  time  t  (assumed  constant 
for  all  t^)  ° 

ETERM  =  exp  (  -T  /  ri  )  (3.17) 

This  concludes  the  development  of  the  dynamics  portion 
of  the  EKF  structure  and  associated  parameter  definitions. 
Later  in  this  chapter  the  method  for  initializing  these 
parameters  with  actual  numerical  values  are  discussed. 

Ill . 4  Measurement  Model 

III . 4 . a  General  Definition.  A  nonlinear  measurement 
model.  Eg. (2.10),  results  from  the  nonlinear  transform 
equations  that  convert  the  noise-corrupted,  orthogonal,  2D 
camera  image  plane  coordinates  to  the  3D  intermediate 
reference  frame  as  detailed  in  Appendix  A. 

For  this  study,  each  camera's  object  plane  size  was  set 
to  a  10.0  x  10.0  ft  square  and  512  x  512  pixel  resolution 
level.  This  numbers  were  selected  purely  because  they  were 
realistic  and  nice  to  work  with.  These  object  plane 

dimensions  result  in  » 0.02  ft  per  side  for  each  pixel  or 

2  .  . 

«0. 0004  ft  per  pixel  m  area.  From  this  resolution 

definition,  the  measurement  noise  is  set  at  a  constant  three- 

sigma  level  of  ±  one  pixel  to  give  a  one-sigma  covariance: 

(t.  )  =  (  0.02  /  3.0  )2  =  4.444  X  10"4  ft2  (3.18) 

nun  iv 
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I 


I 


i 


Since  this  design  involves  two  cameras  that  each  give  a 


2D  position  pair,  (x1,y1)  and  (y.  ,y  )  ,  the  4D  measurement 
vector,  .^(t^) ,  is  now  defined  as  simply  augmenting  the  two 
camera  output  pairs  as: 

*T‘tk>  =  [X1  yi  X2  y2Il  <3'19 


From  Eqs. (A-5)  through  (A-8)  of  Appendix  A,  the 
resulting  4D  vector  of  functions  h[x(t~),tk]  becomes: 


where 


2  D  ft  (  R  x4  +  D  ) 

2  (  *4  x?  +  x24  +  D2  )  +  R  D  (  3  x4  +  Z7  ) 

R  D  (  x?  -  x4  ) 

R  (  *7  +  *4  )  +  2  D 

R  D  H  (  x  +  x  ) 

- - - - - - -  -DX- 

R  (  -  X4  )  +  2  D 


D  -  R  X. 

4 

R  D  (  x?  +  x4  ) 

R  (  “  >^4  )  +  2  D 

(3.20) 


R  =  (  2.0  )1/'2 

D  =  the  camera  focal  length  defined  as  the  distance 
from  the  intermediate  frame  origin  to  each  camera 
lens  (set  to  20  x  21'  per  Section  I.l.b.l). 


This  leaves  the  four-by-nine  matrix  H[x(tk) ,tk]  and  each  of 


its  elements  defined 

by  Eq. (2.12) 

to 

be: 

Hl,l 

0 

0  Hl,4 

0 

0 

Hl,7 

0 

0 

IS 

i 

• 

■ 

LJ 

II 

0 

0 

0  H2,4 

0 

0 

H2 , 7 

0 

0 

H3,l 

0 

0  H3 , 4 

0 

0 

H3 , 7 

0 

0 

0 

0 

0  «4,4 

0 

0 

H4 , 7 

0 

0 

(3.21) 

where,  in  terms  of  the  following  three  denominator  terms: 


DENA  =2  (  *4  X?  +  ft2  +  D2  )  +  R  D  (  3  x4  +  X?  ) 


DENB  =  R  (  ft?  +  X4  )  +  2  D 


(3.22) 

(3.23) 


DENC  =  R  (  ft?  -  x4  )  +  2  D 
the  entries  of  Eg. (3.21)  may  be  evaluated  as: 


(3.24) 


H 


1,1 


2  D  (  R  x4  +  D  ) 
DENA 


(3.25) 


H 


3,1 


D  R  (  X  +  X  ) 

- L - 2 -  _  D 

R  (  X  -  X  )  +  2  D 


D  -  R  X, 


(3.26) 


H 


1,4 


2  D  R  ft. 


DENA 


2  D  k1  (  R  ft4  +  D  ) 
DENA2 
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(2*7+4*4+3DR) 


DENB 


2  (  x?  -  *4  ) 
DENB2 


2  D  X1  (  x?  +  *4  ) 
DENC 


-DR*. 


(  D  -  R  X4  ) 


D  R  xx  2  D  x1  (  x?  +  X4  ) 
DENC  DENC2 


D  -  R  X, 


DR  2  D  (  x?  +  x4  ) 
DENC  DENC2 


“2  D  (  2  S.  +  D  R  )  (  R  St.  +  D  ) 


DENA 


DR  2  D  (  -  *4  ) 


2,7  = 


DENB 


DENB 


(3.27) 


(3.28) 


(3.29) 


(3.30) 


(3.31) 


(3.32) 
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«• 


D  R  ^  2  D  (  X?  +  X4  ) 

DENC  DENC2 


D-R4. 

4 

DR  2  D  (  X?  +  X4  ) 
DENC  DENC2 


(3.33) 


(3.34) 


The  measurement  corruption  noise  strength  matrix,  R(tk) 
has  no  off-diagonal  terms  since  each  measurement  direction  is 
mutually  decoupled  by  the  orthogonal  measurement  geometry. 
Eqs.(2.ll)  and  (3.18)  then  define  R(tk)  as: 

4.44  x  10-4  000 

0  4.44  X  10"4  0  0 

0  0  4.44  X  10-4  0 

000  4.44  X  10“4 

(3.35) 

2 

where  R(tk)  has  units  of  ft  and  is  modelled  as  constant  for 
all  tk. 

Ill . 4 .b  Measurement  Ambiguity. 

Two  CCD  cameras  used  to  collect  separate  orthogonal  2D 
position  measurements  may  seem,  at  first  thought,  to  be 
angle-only  position  measurements  with  no  range  measurement 
information.  However,  that  is  not  the  case.  The  3D  position 
result  from  the  measurement  geometry,  Eqs. (A-2)  through  (A-4) 


B(tk)  = 
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of  Appendix  A,  does  include  the  necessary  range  information 
for  discrimination  between  multiple  fragment  positions.  This 
is  true  as  long  as  the  measurements  of  the  two  cameras  are 
properly  correlated.  Therefore,  the  problem  exists  as  an 
observability  ambiguity  for  multiple  measurement 
combinations,  such  that  the  EKF/JPDA/RTS  package  does  not 
know  "a  priori”  which  fragment  in  camera  #1  correlates  to 
which  fragment  in  camera  #2.  Once  this  ambiguity  is 
overcome,  full  3D  observability  exists.  This  ambiguity  is 
now  illustrated  for  no  clutter,  perfect  detection  situations: 


Let:  fragment^  represent  a  measurement  coordinate  pair  (x,y) . 
Case  1:  Simple  one  fragment  camera  images: 


Figure  3.3.  One  Fragment  Image  Example 


For  the  above  images  from  cameras  #1  and  #2  there  is  no 
confusion  that  fragment.^  and  fragment2  are  the  same  fragment 
(assuming  perfect  detection  and  no  clutter) .  The  resulting 
augmented  measurement  as  defined  by  Eq. (3.19)  becomes: 
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=  [  (fragment^,  (fragment^  ] 


(3.36) 


Case  2 :  Given  the  following  multiple  fragment  camera  images: 


Camera  #1 

Camera  #2 

♦  fragment1 

♦  fragment3 

♦  fragment2 

♦  fragment^ 

Figure  3.4.  Multiple  Fragment  Image  Example 


From  Figure  3.4,  all  that  can  be  determined  immediately  is 
that  there  are  two  fragments.  There  is  no  information 
available,  for  example,  which  of  the  two  measurements  in 
camera  #2  matches  to  fragment1  in  camera  #1.  The  result  is 
two  possible  measurement  combinations: 

z.  =  [  (fragment.),  (fragment  )  ]  (3.37) 

or  J 

z2  =  [  (fragment^,  (fragment^,)  ]  (3.38) 

The  same  problem  occurs  for  the  second  measurement  of  camera 
#1,  giving: 

z_  =  (  (fragment,),  (fragment  )  ]  (3.39) 

or 

Z4  =  [  (fragment2),  (fragment^  ]  (3.40) 

The  final  result  is  that  given  p  camera  #1  measurements  and 

q  camera  #2  measurements,  where  p  =  q,  there  are  p*q 
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candidate  measurements  for  a  perfect  detection/no  clutter 
assumption.  However,  such  an  assumption  is  unrealistic  for 
the  arena  test  environment;  therefore,  p  might  not  equal 
q.  For  such  cases,  p  and  q  may  each  be  greater  than  or 
less  than  the  true  number  of  fragments.  This  further 
increases  the  total  possible  number  of  measurement 
permutations  and  combinations.  For  a  full  description  of 
this  event  algebra,  see  Appendix  B. 

From  the  previous  illustration,  and  the  derivation  of 
Appendix  B,  one  quickly  sees  that  the  number  of  possible 
candidate  fragment  measurements  grows  exponentially  as  the 
number  of  true  or  clutter  fragment  measurements  increase. 

The  selection  of  the  correct  measurement  combinations  from 
the  candidate  set  gives  the  desired  3D  position  information 
for  each  fragment.  Fortunately,  the  JPDA  algorithm 
eliminates  or  reduces  the  effect  of  the  probabilistically 
unlikely  candidate  measurements  primarily  through  the 
validation  gate  process  (as  in  Figure  2.8)  and  secondarily 
though  the  joint  probability  residual  weighting.  Figure  3.5 
gives  a  general  system  block  and  signal-flow  diagram  of  the 
arena  test  EKF/JPDA/RTS  processing  system. 
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Figure  3.5.  Full  System  Block  and  Signal-Flow  Diagram. 


Ill . 5  Filter  Initialization  Procedure 


The  transient  performance  of  the  EKF  during  acquisition 

is  largely  determined  by  the  initial  estimates  of  state 

values,  x(tQ) ;  initial  covariance  matrix  values,  P(tQ) ; 

dynamics  driving  noise,  Qd(to);  and  the  dynamics  model 

parameters  for  the  acceleration  time  constants  (r  ,  r  and 

x  y 

r  ) .  The  closer  these  initial  values  are  to  the  true  values, 
the  smaller  the  transient  and  associated  tracking  errors. 

The  following  sections  present  two  possible  methods  to  obtain 
the  above  initialization  values  and  parameters,  1)  by 
statistical  data  method,  and  2)  by  exploiting  "a  priori" 
arena  geometry  information.  The  latter  of  the  two  methods 


being  the  choice  for  implementation. 
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Ill • 5 . a  Statistical  Initialization  Approach.  The 
fragment  cloud  that  disperses  from  a  detonated  munition  can 
be  reasonably  assumed  to  include  fragments  of  various  sizes 
(cross  sectional  area)  and  various  speeds  (velocity 
magnitude) .  One  approach  for  selecting  initial  values  is  to 
compute  the  statistically  averaged  values  for  velocities  and 
acceleration  time  constants  (ATCs)  from  previous  test  data 
and  selecting  abnormally  high  covariance  and  dynamics  noise 
values.  The  high  E(tQ)  and  Qd(tQ)  values  initially 
increase  the  Kalman  gain  such  that  the  increased  filter 
bandwidth  quickly  (hopefully)  dampens  out  any  transients. 

However,  such  a  method  requires  that  previous  test  data 
exist.  In  addition,  because  of  the  large  range  of  velocities 
and  acceleration  time  constants  that  may  occur  in  one  test- 
blast,  many  of  the  acquired  fragments  will  have  true  state 
values  far  from  the  expected  averages.  This  results  in  large 
start-up  transients  or  worse — loss  of  track. 

For  this  research  effort,  an  alternate  initialization 
method  is  developed  and  gives  excellent  results.  The 
deterministic  arena  set-up  geometry  with  its  known  dimensions 
provides  significant  "a  priori”  information  that  the 
initialization  process  may  exploit. 

Ill . 5 .b  "A  Priori”  Information.  The  deterministic 
geometry  model  may  be  exploited  by  one  obvious  assumption: 
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Since  the  test  munition  is  detonated  at  the  W-frame 
origin,  there  is  a  high  probability  that  any 
fragment  that  enters  a  camera’s  field-of-view  (FOV) 
initiated  its  trajectory  from  the  W-frame  origin  at  time 

V 

Recall  in  Chapter  II  that  first  light  sensors  detect  the 

breakup  of  the  case  of  the  test  munition  to  start  the  system 

clock  at  tQ.  As  a  fragment  enters  a  camera’s  FOV,  the 

fragment's  measured  diameter,  position  and  the  time  elapsed 

since  t  can  be  used  to  calculate  that  fragment's  initial 

position,  velocity,  drag  induced  acceleration,  and 

acceleration  time  constants.  From  these  values,  r  ,  r  ,  r  , 

x  y  z 

x(to) ,  P(tQ) ,  and  Q(tQ)  may  be  generated  to  initialize  the 

EKFs .  This  process  will  now  be  developed. 

The  initialization  procedure  begins  by  transforming  the 

initial  measured  position,  z(tm) ,  into  the  W-frame  (Appendix 

A)  to  become  E(tm) •  The  average  speed  over  the  time 

interval,  v  „  rt  ,t  ],  is  simply  the  magnitude  of  the  W- 
avg1-  o  m 

frame  position  vector,  I  >  divided  by  the  elapsed  time, 

t  -t  . 
m  o 

Begin  by  setting: 

E( tQ)  =  0  (at  W-frame  origin) 

t  *  0  (time  of  detonation) 

t  =  the  time  of  first  measurement 
m 

where 

l»<Vl  -  t  Px<V  +  Py<V  +  Pz<V  51/2  <3-41> 

and 


I 


I 


i 


i 
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To  avoid  requiring  three  separate  sets  of  calculations  for 
each  axis  direction,  the  following  procedures  are  carried  out 
once  with  the  scalar  magnitudes  from  above.  The  obtained 
scalar  results  will  then  be  projected  back  into  the 
appropriate  axis  components.  From  Eq.(3.4),  the 
deterministic  portion  of  the  acceleration  model,  in  scalar 
form,  may  be  extracted  as: 

a (t)  »  -l/r  a (t)  (3.43) 

and  taking  the  time-derivative  of  Eq.(3.3b)  for  a  scalar  case 
where  |v(t) |  is  replaced  with  v(t)  gives: 


aD(t)  = 


2  K1  V  ( t )  V  ( t ) 


(3.44) 


where 


v (t )  =  aD(t) 


(3.45) 


Substituting  Eqs.(3.45)  and  (3.44),  along  with  vaVg[t0'"tjn] 
and  d(tm)  in  place  of  v(t)  and  d  respectively,  into 
Eq.(3.43)  and  solving  for  r  as  r  gives: 


d<v 


avg 


2  K1  w  w 


(3.46) 


Note  in  Eq.(3.46),  that  TaVg  is  dependent  upon  the  measured 

fragment  diameter,  d(t  ),  and  v _ [t  ,t  ].  Under  the 

J  m  avg  o  m 


assumptions  that  TaVg  for  a  fragment  remains  near  constant 
and  that  gravitational  acceleration  is  negligible  for  all 
time  of  interest,  this  ravg  parameter  is  now  placed  into  a 
scalar  continuous-time  form. 


Note  that  vector  notation  will  be  temporarily  dropped,  where 
now  indicated  position,  velocity,  and  acceleration  values  are 
scalar  magnitudes  of  the  previous  vectors: 


p(t)  =  aD(to)  [  ~TaVg  +  t  Tavg  +  ravg  exp(-t/Tavg)  ] 

+  t  v(to)  +  p(to)  (3.47) 


v(t)  =  aD(to)  ravg  [  1  -  exp(-t/ravg)  ]  +  v(tQ) 

aD(t)  =  W  exP<"t/ravg) 

The  above  equations  are  all  dependant  upon  aD(tQ) • 
may  be  found  through  Eq.(3.3b)  as: 


(3.48) 

(3.49) 
which 


aD(tO>  = 


K1  V<V 


(3.50) 


In  this  equation,  fragment  diameter  is  constant  so  d(tm)  = 

d(t  )  =  d,  and  v(t  )  is  the  solution  of: 
o  o 

nM.  ,  _  K1  [  fcm  Tavq  “  ravo  +  ravq  exP^VTavq)  1  v(to)2  . 

P(  m)  =  d 

+  tm  v(tQ)  (3.51) 

where 


A  = 


[  ^m  Tavg  Tavq  +  Tavq  eX^  Si/Tavq^ 


(3.52) 
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B  =  t 


c  -  -P(V 


(3.53) 

(3.54) 


and  now  solving  the  quadratic  below  gives  v(tQ): 


v(to)  - 


-B  +  (  B2  -  4  A  C  )1/2 


(3.55) 


Note  that  only  the  positive  root  solution  for  Eq. (3.55) 
is  applicable  since  scalar  magnitudes  are  involved. 

Once  v(tQ)  and  aD(tQ)  have  been  obtained  by  solving 
Eq.(3.55)  and  substituting  back  into  Eq.(3.50),  v(tm)  and 
aD(tm)  may  be  found  through  Eqs.(3.48)  and  (3.49)  evaluated 

at  V 

v<V  '  W  f avg  '  1  ‘  ex?<  ‘S.  '  'avg  >  >  +  v<V  (3'56> 


W  *  W  exP<  -‘m  '  Tavg  > 


(3.57) 


The  scalar  magnitude  values  of  raVg,  v(tm),  and  a0(tItl) 
are  then  converted  into  appropriate  X w,  Yw,  and  Zw 
component  vectors  by  transformation  terms,  XTRM,  YTRM,  and 
ZTRM  solved  by  trigonometric  projections  to  give 

Tx  -  XTRM  r  ;  vx(tm)  =  XTRM  v(t„)  ;  =  XTRM  «D<V 


rv  '  yTRM  Tavg  vy<V  =  YTR"  V(V  '  aDy(tm>  "  iTRM  W 


'a  =  ZTRM  »  ;  vz(tn)  -  ZTRM  v<tB)  ;  aDz(tm)  =  ZTRM  aD(tm) 


where  the  projection  terms  onto  the  X^  Y^,  and  Zw  axes  are: 
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(3.58) 


XTRM  =  COS (ang1)  SIN(ang2) 


YTRM  =  SIN(ang1)  SIN(ang2) 

(3.59) 

ZTRM  =  pz(tm)  /  |E(tm)| 

(3.60) 

-1 

angi  =  TAN  “(py(tn)  /  Px(tB)) 

(3.61) 

ang2  =  COS  tPz(tm)  /  lB(tn)|] 

(3.62) 

Up  to  this  point,  all  accelerations  have  been  due  to  drag,  so 

now  a  gravitational  acceleration  component  is  included  by 
.  2 

summing  -32.1741  ft/sec  to  the  W-frame  acceleration 
vector's  V?z  direction.  The  final  set  of  required 
transformations  takes  the  derived  initial  conditions  in  the 
W-frame  and  converts  them  to  the  I-frame  (Appendix  A)  where 
the  EKF  operates. 

To  complete  the  EKF  initialization  process,  the  above 

calculated  a_  (t  ),  a_  (t  ) ,  and  a_,(t  )  values  are  squared 
Dx  m'  Dy  m  Dz  m  ^ 

.  .  2  2 
and  inserted  as  the  mean-squared  accelerations  cr  .  a  ,  and 

ax  ay 

2 

a  into  Eqs.(3.9)  through  (3.16)  to  generate  P(t  )  and 
fld(tQ) »  where  Qd(*)  remains  constant  with  time.  Initial 
EKF  test  runs  indicate  not  to  rely  directly  on  the  values  of 

aDx^tm)  '  aDy(tm)'  and  aDz(tm)  to  set  the  W*  This  is 
because  the  above  acceleration  drag  terms,  aD(*),  may 

calculate  close  to  zero:  resulting  in  near-zero  Q^(') 

values.  This  allows  the  EKF's  gain  matrix  to  decrease  too 

rapidly  and  possibly  breaking  track  from  the  true  trajectory. 
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The  solution  to  this  initialization  shortfall  is  to 


incorporate  a  minimum-level  drag  acceleration  value,  aDxnin' 
that  was  set  (by  trial-and-error  tuning)  to  one  third  the 
value  of  aDX(tjn)  since  the  X^.  direction  is  expected  to 
have  the  highest  aD(*)  component.  This  minimum-level  value 
is  used  to  calculate  fld(tQ)  and  P(tQ)  as  a  default  value 
should  the  above  calculated  aD(*)  values  fall  below  the 

a_  .  threshold. 

Dmin 

This  completes  the  entire  EKF  initialization  process  for 
*(tQ) ,  P(tQ) ,  and  fl(tQ)  required  by  the  Singer  approach. 

Ill . 6  Summary 

This  chapter  has  introduced  the  applicable  system 
reference  frames,  and  given  detailed  equation  structures  for 
the  applied  truth  model,  linear  EKF  Singer  dynamics  model, 
and  nonlinear  EKF  measurement  model.  In  addition,  a 
detailed  break-down  of  the  applied  EKF  initialization 
procedure,  based  on  known  arena  geometry  and  the  initial 
fragment  measurement,  was  given. 


! 
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IV.  SIMULATION 

This  chapter  presents  the  structure  and  implementation 
of  the  fragment  tracking  digital  computer  simulation 
developed  in  this  study.  The  simulation  was  written  entirely 
in  double  precision  FORTRAN  77  and  is  self-supporting  (i.e. 
no  IMSL  or  library  system  calls)  on  most  large  mainframes. 

The  sole  exception  is  the  use  of  MATRIX^  (32)  to  generate 
plots.  Overall,  the  simulation  is  divided  into  two 
independent  main  programs,  1)  a  fragment  trajectory  generator 
(truth  model)  and  2)  a  self-contained  EKF/JPDA/RTS 
tracker/ smoother  program.  The  trajectory  generator  program 
creates  data-files  for  the  tracking  program  to  use  as 
"measurements"  for  its  operation. 

IV. 1  Fragment  Trajectory  Generator 

The  fragment  trajectory  generator,  is  a  fully  self- 
contained  program  that  calculates  the  time-varying 
acceleration,  velocity  and  position  of  up  to  six  fragments 
simultaneously  from  a  set  of  initial  conditions.  Four 
"trackable"  fragments,  numbered  one  through  four,  and  two 
"clutter"  fragments,  numbered  five  and  six  may  be 
selectively  generated.  The  resulting  data  is  stored  in  a 
formatted  datafile  to  provide  "measurements"  for  subsequent 
processing  by  the  tracker/ smoother  program  as  well  as  an  easy 
to  read  "proof"  file  that  lists  each  fragment's  time  varying 
position,  velocity,  and  acceleration  state  values  in  the  W- 
frame,  I-frame,  and  C-frame. 


The  fragment  trajectory  dynamics  are  based  on  the 
fragment  acceleration  truth  model,  Eqs.(3.1)  through  (3.3b) 
presented  in  Chapter  III.  Fragment  trajectories  are  modelled 
in  accordance  with  Appendix  A  and  are  initiated  from  the 
arena  center  (W-frame  origin) ,  fi(tQ)  =  0,  and  an  initial 
velocity  vector,  v(tQ) .  The  acceleration,  a[v(t),g],  on  each 
fragment  is  defined  as  the  velocity  magnitude  (speed) 
dependent  atmospheric  drag  component  (initialized  with  v(tQ)) 
summed  with  a  constant  gravitational  component,  g.  Using  the 
two  integration  steps,  as  illustrated  in  Figure  3.2,  each 
fragment's  acceleration  vector,  is  recursively  integrated 
over  one  frame  period  (200  jisec)  to  produce  the  velocity 
vector,  v(k+l) .  This  new  velocity  is  used  to  both  calculate 
the  new  acceleration,  &[X(K+l),g],  and  is  integrated  to 
produce  the  position  vector,  p(k+l) . 

Integration  is  accomplished  by  the  trapezoidal  technique 
where  each  integration  period  (one  frame  period)  is  divided 
into  1000  trapezoidal  summations.  The  1000  summations 
quantity  was  selected  assuming  the  near  straight-line 
trajectory  would  be  adequately  smooth  for  this  resolution. 

Constant  values  programmed  into  the  fragment  generator 
program  include: 


Atmospheric  density: 
Fragment  density: 
Fragment  diameter:  d 
Coefficient  of  drag:  C_ 


0.076474  Ibm  /  ft  (sea  level) 
491.0  lbm  /  ft3  (low-carbon  steel) 
0.5  inches  (0.0416667  ft) 

1.0  (as  assumed  in  Chapter  1) 
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Camera  to  W-frame  origin  vertical  displacement:  5.0  ft 


The  following  default  values  may  be  interactively  changed  by 
the  user: 

Arena  center  to  camera  location  radius:  Rd  =  50.0  ft 

Camera  focal  length:  D  =  20* (2. O)1^2  =  28.28427712474  ft 

I-frame  origin  angular  displacement  from  X w:  0  =  30.0° 

Camera  angular  f ield-of-view  (FOV) :  a  =  20.0499757242° 

(This  gave  the  "nice”  10.0  ft  image  plane  size  below.) 

The  default  parameters  result  in  the  following  dimensions: 

W-frame  origin  to  I-frame  origin  vertical  displacement: 

Sz  =  15.0  ft 

Camera  object  plane  size:  Size  =  10.0  ft 

The  initial  fragment  speed  (velocity  magnitude)  is 
interactively  prompted  and  entered  by  the  user  when  the 
program  is  started.  Each  fragment’s  initial  direction  is 
pre-determined  by  a  set  of  polar-coordinate  angles  such  that 
the  trackable  fragments,  one  through  four,  pass  through  the 
cameras’  FOV  intersection  (coverage  volume).  Clutter 
fragment  number  five  passes  through  camera  #2’s  ”near-field” 
FOV  while  clutter  fragment  number  six  passes  through  camera 
#l's  ”far-field"  FOV.  These  result  in  fragment  five  being 
only  visible  to  camera  #2  and  fragment  six  only  visible  to 
camera  #1;  creating  an  extra  measurement  in  each  camera  that 
cannot  be  tracked;  or  clutter. 

An  additional  user  selectable  feature  is  the  option  to 

6' 


add  "curving"  accelerations  to  trackable  fragments  two, 
three,  and  four.  These  curving  accelerations  cause  fragments 
two  and  three  to  have  effectively  "crossing  tracks,"  while 
the  curving  acceleration  applied  to  fragment  four  causes  it 
to  curve  downward  sharply  (-Zw  and  -Zj  directions) . 

The  multitude  of  fragment  trajectory  combinations 
possible  for  this  program  to  simulate  allows  thorough  testing 
of  the  EKF/JPDA/RTS  package  from  "easy"  single-fragment 
straight  trajectories,  to  relatively  "difficult"  multi¬ 
fragment,  multi-clutter,  crossing  tracks,  and  curving 
trajectory  scenarios. 

IV. 2  EKF/JPDA/RTS  Tracker/Smoother  Program 

The  EKF/JPDA/RTS  program  reads  the  datafile  generated  by 
the  fragment  generator  program  exactly  as  though  the  file 
were  pre-processed  image  measurements.  These  measurements 
include  for  each  camera:  the  number  of  fragments  present  in 
the  image,  the  centroid  (x,y)  coordinates  for  each  observed 
fragment,  and  the  diameter  of  each  observed  fragment.  These 
truth  model  measurements  are  provided  in  both  the  four¬ 
dimensional  measurement  space  as  zft.^)  and  as  three- 
dimensional  position,  velocity,  and  acceleration  state 
values,  H(t^) ,  in  the  I-frame  for  subsequent  error  and 
statistical  analysis  of  the  algorithm’s  performance. 

IV. 2 . a  Measurement  Corruption  Noise  Technique.  The 
truth  model  measurements  (fragment  position  and  diameter)  are 
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corrupted  by  the  summation  of  a  pseudo-random  number  noise 
sequence  to  the  uncorrupted  measurement  vector  (see  Figure 
4.1).  The  unity  strength  noise  generator  shown  in  Figure  4.1 
includes  a  pseudo-random  number  generator  from  which  twelve 
independent  outputs  are  summed,  followed  by  subtracting  six 
from  this  accumulated  sum  for  each  element  of  a  four¬ 
dimensional  noise  vector.  This  results  in  a  6D  (4D  position 
+  2  diameter  measurements)  approximate  zero-mean  Gaussian 
distributed  noise  vector  of  covariance  equal  to  a  six-by-six 
identity  matrix.  The  overall  output  variance  of  this  noise 
vector  is  adjusted  by  multiplying  each  covariance  term  by  the 

desired  standard  deviation,  a  ,  then  summed  to  the  4D  truth 

mm 


Figure  4.1.  Measurement  Noise  Corruption  Process 


measurement  vector,  ^ciean'  and  fra9ment  diameters,  dcan,i 
and  dcam2  to  give  a  noise-corrupted  measurement  vector  and 
diameter  measurements  for  processing  by  the  EKF/JPDA/RTS 


tracker/ smoother  package. 
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IV. 2 .b  Tracker /Smoother  Specifics.  This  program  allows 
the  user  to  adjust  the  following  parameters  at  the  beginning 
of  each  simulation  run: 

The  truth  model  fragment  diameter  measurement  noise  variance. 

The  truth  model  fragment  position  measurement  noise  variance. 

The  JPDA  algorithm  validation  gate  size:  g,  of  Eg. (2.18). 

The  JPDA  algorithm  fragment  probability  of  detection:  Pn,  of 
Eq.(2.20)  (all  fragments  are  assumed  to  have  equal  UPD) . 

The  JPDA  algorithm  Poisson  clutter  density:  C,  of  Eq.(2.20). 

For  all  runs  presented  in  this  study,  both  of  the  above 

truth  model  noise  levels  where  set  equivalent  to  the  filter's 

expected  measurement  noise  levels  of  one  pixel  or 
-4  2 

4.444  x  10  ft  .  This  noise  level  was  selected  simply  as  an 

apparently  logical  place  to  start  any  analysis.  The 

selection  of  the  validation  gate  size  was  rather  arbitrary, 

using  trial  and  error  to  arrive  at  reasonable  values.  Values 

of  g  =  3.0  gave  good  tracks  for  straight  fragment 

trajectories  while  g  =  4.0  became  necessary  to  maintain 

reliable  tracking  for  curving  trajectories.  Since  the 

validation  gating  process  takes  place  in  the  4D  measurement 

4 

space  of  units  ft  ,  the  volume  represented  by  g  is  a  4D 

"statistical"  volume  that  cannot  be  easily  interpreted  into 

.  3 

3D  real-world  units  of  ft  . 

The  probability  of  fragment  detection,  PD,  was 
arbitrarily  set  at  0.95,  and  the  Poisson  clutter  density,  C, 
was  set  to  unity  for  all  presented  simulation  runs.  These 
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values  were  chosen  as  "nice"  values  work  with  and  due  to 
lack  of  fragment  data  to  select  otherwise. 

Similar  to  the  fragment  generation  program,  this  program 
contains  the  arena  dimension  and  camera  location  constants 
necessary  to  define  the  measurement  geometry  and  perform  the 
required  transformations  between  any  of  the  reference 
frames.  These  dimensions  are  especially  important  in  the  EKF 
initialization  process,  since  this  process  relies  heavily  on 
arena  dimensional  constants  in  its  calculations. 

Once  all  necessary  optional  inputs  are  supplied  to  the 
tracker/ smoother  program,  it  processes  the  data  as  outlined 
in  Chapters  II  and  III.  The  resulting  fragment  state 
estimates  for  both  the  EKF  and  the  RTS  smoother  outputs  are 
compared  to  the  uncorrupted  truth  model  states  to  generate 
the  time-varying  estimation  errors,  x(*)  -  x(*)»*  filter 
computed  standard  deviation  values,  of;  and  true  error 
standard  deviation  values,  at  ;  for  Posit;i-on'  velocity,  and 
acceleration  in  the  Xj,  Yj,  and  ZI  directions.  These 
statistical  values  are  then  formatted  and  written  out  into 
two  separate  MATRIX^  loadable  datafiles  for  plotting.  The 
first  datafile  contains  the  EKF  performance  data  and  the 
second  contains  the  RTS  smoother  performance  data.  A 
representative  example  plot  is  shown  in  Figure  4.2.  The 
solid  center  line  is  the  estimate  error,  e;  the  top  solid  and 
bottom  dotted  lines  that  are  symmetric  about  zero  are  the 
filter's  standard  deviation,  ±  af;  and  the  top  dotted  line 
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and  bottom  dashed  line  symmetric  about  the  estimate  error 
line  describe  the  true  one-sigma  standard  deviation  bound  of 
the  estimate  error,  e  ±  ;  where  each  of  these  values  are 

calculated  in  the  vector  sense  by: 


e(-) 


ef  ( * ) 


1 
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1 
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£(•) 


diag 


n 


true  states 
estimate  states 

diagonal  elements  of  the  filter 
estimate  covariance  matrix 

number  of  Monte-Carlo  runs 


(4.1) 


(4.2) 


(4.3) 


IV. 3  Summary 

This  chapter  has  presented  the  characteristic 
descriptions  of  the  FORTRAN  fragment  generation  and 
tracker/ smoother  simulation  programs  developed  in  this 
study.  This  was  followed  by  a  detailed  description  of  the 
system's  applied  noise  corruption  technique.  The  chapter  was 
completed  with  a  description  of  the  simulation  generated 
information  displayed  in  its  output  plots;  including  an 
example  plot. 
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m  dx 


Figure  4.2.  Example  Performance  Plot 


V.  ANALYSIS 


VI . 1  Simulation  Runs 

Three  major  trajectory  categories  were  run  against  the 
presented  EKF/JPDA/RTS  package  by  Monte-Carlo  simulation: 

Category  l.)  Single  fragment:  straight  trajectory  with 
no  clutter. 

Category  2.)  Multiple  simultaneous  fragments  (three): 

each  with  straight  trajectory  and  no 
clutter. 

Category  3.)  Multiple  simultaneous  fragments  (five): 

one  with  straight  trajectory,  two  with 
curving  trajectories  that  cross,  and  two 
extraneous  clutter  fragments. 

Each  of  the  above  trajectory  categories  were  run  at 
three  initial  fragment  speeds,  |v(to) | ,  of  3000,  6000,  and 
10000  ft/sec  (total  of  nine  simulation  runs) .  These  runs 
produced  performance  plots  for  each  modelled  fragment’s  EKF 
and  RTS  smoother  tracking  errors  in  position,  velocity,  and 
acceleration  for  each  of  the  Xj,  Yj,  and  Zj  directions  (18 
plots  per  fragment) .  A  total  of  324  simulation  performance 
plots  where  generated  from  these  runs  and  are  displayed  in 
Appendix  C. 

The  desired  number  of  Monte-Carlo  runs  for  each 
simulation  was  initially  set  to  100;  howeve-r,  computational 
run-time  constraints  required  this  number  be  reduced  to  20 
runs  per  simulation.  This  degraded  the  overall  "smoothness" 
of  the  plots  but  still  allowed  the  important  performance 
trends  to  stand  out.  Example  plots  of  20  versus  100  Monte- 
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Carlo  run  simulation  comparisons  are  shown  in  Figures  5.1  for 
the  EKF  output  and  Figure  5.2  for  the  RTS  smoother  output. 


Figure  5.1.  EKF  Plots  for  20  and  100  Monte-Carlo  Runs. 
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The  previous  example  plots  serve  to  emphasize  that  the 
presented  simulation  results  (in  plots  and  tables)  would,  in 
general,  show  better  performance  if  the  100  runs  could  have 
been  applied  in  all  cases. 


V. 2  Performance  Summaries 

V. 2 .a  Category  One.  Position  and  velocity  tracking 
performances  were  excellent.  This  was  attributed  largely  to 
the  "expected  straight  path”  assumption  imposed  by  the  filter 
initialization  routine  and  the  single  fragment,  with  its  sole 
measurement,  allowing  the  JPDA  algorithm  to  remain  dormant. 
The  obtained  position  estimation  accuracies,  even  with  a 
three-sigma  measurement  noise  level  of  ±  one  pixel,  often 
exceeded  the  one  pixel  resolution  level  (K0.02  ft)  of  the  CCD 
cameras  by  an  order  of  magnitude.  This  "better  than 
expected”  performance  raised  suspicion  that  a  coding  error 
may  exist  in  the  pseudo-random  noise  generator,  causing  lower 
measurement  noise  corruption  levels  than  desired. 

Therefore,  additional  coding  was  added  to  calculate  and  plot 
the  true  standard  deviations  for  each  of  the  following: 


1. )  Each  noise  corrupted  element  of  the  4D  position 

measurement  vector,  z(*). 

2. )  Each  camera's  noise  corrupted  diameter  measurement, 

dcaml  and  dcam2  * 

3. )  The  equivalent  3D  noisy  position  components  of  the 

state  vector,  x(*),  resulting  from  transforming 
the  4D  s(*)  vector  to  the  x(’)  space. 


An  additional  simulation  run  was  made  for  a  3000 
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ft/sec,  category  one  fragment,  with  the  position  and 
measurement  noise  levels  set  for  a  three-sigma  standard 
deviation  level  of  one-pixel  (  0.02  ft  )  as  it  was  in  all 
previous  simulation  runs.  This  defined  the  expected  one- 
sigma  value  to  be  0.02  /  3  k  0.0067  ft.  In  addition,  1000 


Monte-Carlo  runs  were  chosen  to  ensure  statistically  steady- 
state  results  due  to  the  ’000  •'sample"  size.  The  top  plot 
of  Figure  5.3  shows  the  resulting  standard  deviations  for 
each  element  of  the  noise  corrupted  z(*)  vector,  while  the 
bottom  plot  shows  that  for  each  camera's  noise  corrupted 
diameter  measurement.  The  plotted  levels  of  a  »  0.0067  ft 
validity  the  measurement  noise  corruption  process  and  further 
validity  the  category  1  and  2  results.  However,  an 
interesting  system  characteristic  was  discovered  when 
plotting  the  standard  deviations  of  the  equivalent  noises 
into  the  x(*)  space,  as  shown  in  Figure  5.4.  This  plot 
indicates  that  even  though  each  direction  of  the  4D  jz(*) 
vector  have  equal  noise  levels  (as  in  Figure  5.3),  the 
nonlinear  transformation  into  the  3D  position  components  of 


x(*)  results  in  unequal  noise  levels  for  the  Xj,  ,  and  Zj 
directions.  This  plot  indicates  the  Xj  direction  having 
the  least  noise  sensitivity  while  the  Yj  direction  having 
the  most  (worse)  noise  sensitivity. 
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STANDARD  DEVIATION  (FT)  STANDARD  DEVIATION  (FT) 


Figure  5.3.  True  Standard  Deviation  Noise  Level  Plots  from 
Desired  Level  of  «  0.0067  ft  and  1000  Monte- 
Carlo  Runs. 
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Figure  5.4.  Resulting  Noise  Level  Standard  Deviations  for 

the  X  ,  Y  ,  and  Zj  Axis  Directions  from  Noise 
Levels'1  Indicated  in  Figure  5.3. 


V.2.a. 1  Convergence.  The  category  one  filter 
estimate  covariances,  P(*),  for  position  and  velocity, 
converged  quickly  and  correlated  well  with  the  true  error 
covariances  to  indicate  reasonably  good  "tuning"  of  the  EKF’s 
dynamics  noise  strength,  Qd(*)* 

In  contrast,  the  filter  acceleration  state  covariances 
(plotted  as  standard  deviations)  in  all  simulation  runs 
demonstrated  poor  correlation  to  the  true  error  standard 
deviations.  In  all  runs,  the  filter  acceleration  covariance 
was  substantially  larger  than  the  true  covariance.  This  was 
mainly  attributed  to  the  long  transient  response  of  the  EKF's 
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acceleration  model.  Typical  acceleration  time  constants,  r, 
were  observed  to  range  from  0.01  to  0.5  second  while  a  10,000 
ft/sec  track  interval  would  only  last  80  nsec.  Therefore, 
the  comparatively  "slow"  acceleration  model  could  not  be 
converged  quickly  enough  by  the  "brief"  interval  of 
measurements . 

Table  5.1  summarizes  the  performance  results  taken  from 
the  category  one  RTS  smoother  outputs.  Error  data  presented 
in  this  table  (and  following  tables)  includes  the  worst-case 
error  value  for  the  entire  track  duration,  [tQ,tf]  and  the 
percentage  of  that  error  to  the  true  state  value  at  its  time 
of  occurrence.  For  example,  the  X-axis  velocity  estimate 
error  value  in  Table  5.1  for  a  3000  ft/sec  fragment  was 
determined  by  observing  Figure  C.4  in  Appendix  C  and  noting 
that  the  largest  error  (worst-point)  was  approximately  2.0 
ft/sec  and  occurred  at  time  0.0212  sec.  The  error 
percentage  is  then  calculated  by  dividing  the  above  error 
value,  2.0,  by  the  true  state  value  at  the  worst-point  time 
of  0.0212  sec  (obtained  from  truth  model  computer  listings) 
as  2407.354  ft/sec: 

•  100.0  «  0.083%  (6.1) 

For  cases  where  the  truth  model  value  is  zero  or  close 
to  zero,  Eq.(6.1)  becomes  undefined  or  grows  misleadingly 


large.  An  infinity,  ®,  symbol  is  placed  in  the  table's 
"percent"  square  to  indicate  such  occurrences. 

Regarding  position  estimation  errors,  since  the  true 
position  values  vary  across  an  interval  of  «[-5.0,+5.0]  ft, 
then  expressing  such  errors  as  percentages  would  also  be 
misleading.  Therefore,  percentage  values  are  considered  not 
applicable  for  position  and  an  N/A  is  placed  in  the  tables 
appropriate  "percent"  squares.  For  all  figures  and  tables, 
position  error  values  are  in  feet,  velocity  errors  are  in 
ft/sec,  and  acceleration  errors  are  in  ft/sec  . 

Pointing  out  some  highlights  from  Table  5.1  and  Figures 
C.l  through  C.54  in  Appendix  C  include: 


1. )  No  smoothed  position  estimate  errors  greater  than 

0.002  ft  (0.024  inches). 

2 .  )  No  smoothed  velocity  estimate  errors  greater  than 

0.253%  of  the  true  value. 

3 .  )  Smoothed  acceleration  estimate  errors  vary  across  a 

rather  large  range  from  0.503%  to  24.14%  (a  good 
indicator  of  potentially  unreliable  performance) . 

4. )  The  EKF  initialization  method  devised  in  this 

research  demonstrated  low  sensitivity  to  the 
fragment  diameter  measurement  noise  (critical  to 
derive  r  parameter)  and  performed  well  at  all 
three  fragment  speeds. 


Overall,  the  EKF/JPDA/RTS  package  performed  very  well 
(excluding  acceleration  performance)  for  category  one 
simulation  runs. 
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Table  V.I.  Category  1  (single  fragment)  Worse-Case  Error 
Performances 


For  initial  speed: 

ls<V 

=  3000. 

)  ft/sec 

X-Axis 

- - 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0012 

N/A 

-0.0018 

N/A 

-0.0011 

N/A 

Velocity 

+2.0 

0.083 

-0.085 

00 

.  . 

+  0.50 

0.059 

Accel . 

+1125 

6.45 

+  30.0 

00 

.  

-500 

8.41 

For  initial  speed: 

|v(to) 

=  6000. 

)  ft/sec 

X-Axis 

Y-Axis 

_ 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0008 

N/A 

-0.0019 

N/A 

+0.0010 

N/A 

Velocity 

-1.6 

0.032 

-0.9 

00 

+  2.2 

0.133 

Accel . 

+2000 

2.89 

+800 

00 

+  3400 

13.54 

For  initial  speed: 

lY(t0) 

=  10000 

0  ft/sec 

X-Axis 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.002 

N/A 

— 

-0.0007 

N/A 

-0.002 

N/A 

Velocity 

+7.0 

0.085 

+  0.8 

00 

+7.2 

0.253 

Accel . 

-10150 

0.503 

-125 

00 

-16000 

24.14 

Regarding  overall  steady-state  convergence,  the  small 
number  of  propagate/update  cycles  that  occur  for  tracking 
6000  ft/sec  (twelve  updates)  and  10000  ft/sec  (six  updates) 
fragments  do  not  allow  the  EKF  states  to  fully  converge  down 
to  a  steady-state  condition.  This  becomes  obvious  when 
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comparing  a  3000  ft/sec  run  to  a  10000  ft/sec  run.  The  3000 
ft/sec  run  receives  enough  updates  to  converge  and  "level 
off"  to  a  steady-state  condition  while  the  10000  ft/sec  run 
does  not  receive  enough  updates  to  do  the  same. 

It  was  noted  that  the  EKF's  first-order  Gauss-Markov 
acceleration  model  as  defined  in  Eq.(3.4),  fit  very  close  to 
the  acceleration  truth  model  in  Eq.(3.1)  as  long  as  the 
acceleration  time  constant,  t,  was  correctly  matched  as  a 
function  of  velocity.  This  close  model  matching  becomes 
beneficial  for  cases  of  "no  valid  measurements"  in  the  JPDA 
algorithm  update  cycle.  Such  a  situation  was  observed 
frequently  during  simulation  runs  when  one  or  more  noise- 
corrupted  position  measurements  fell  outside  the  validation 
gate.  This  required  the  EKF  acceleration  model  to  continue 
propagating  without  a  measurement  update  for  one  or  more 
recursion  cycles.  Therefore,  the  closer  the  model  match,  the 
smaller  the  tracking  error  growth  rate  during  an  extended  "no 
valid  measurements"  propagation  period.  Although  the  Singer 
dynamic  acceleration  model  demonstrated  good  matching  to  the 
truth  model  characteristics,  the  overall  acceleration 
estimate  errors  were  considerably  high.  This  was  attributed 
to  the  fact  that  the  acceleration  states  are  two 
differentiations  away  from  the  incoming  noisy  position 
measurements.  The  two  differentiation  steps  accentuate  the 
high  frequency  system  noises  to  levels  that  dominate  the 
acceleration  state  estimates.  This  poor  estimation 
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performance  combined  with  the  previously  discussed  poor 
covariance  performance  leaves  no  choice  but  to  regard  the 
acceleration  estimates  as  unreliable  and  unusable. 

V. 2 . b  Category  Two.  These  multi-fragment  (multi¬ 
measurement)  simulation  runs  are  where  the  JDPA  algorithm  was 
first  put  to  work.  For  these  runs,  the  three  fragments,  with 
their  straight  trajectories,  did  not  pose  any  difficulty 
whatsoever  to  the  JPDA  algorithm,  as  long  as  the  fragment 
trajectories  were  not  too  close  together.  Close  comparison 
of  Figures  C.7  and  C.73  in  Appendix  C  reveals  the  filter 
standard  deviation  lines  in  Figure  C.73  begin  to  diverge 
after  0.0195  seconds.  This  peculiarity  imposed  the  need  for 
a  closer  look  by  running  two  additional  special  runs:  one 
run  for  two  relatively  "close"  fragment  trajectories  (»1.8  ft 
apart)  and  another  run  for  two  relatively  "distant" 
trajectories  (a2.6  ft  apart).  These  scenarios  were  run  for 
100  Monte-Carlo  runs  to  rule  out  if  Figure  C.73  was  only  the 
result  of  its  relatively  small  20  Monte-Carlo  run  sample 
size.  Figure  5.5  displays  two  indicative  plots  that 
summarize  the  results.  The  top  plot  shows  the  X-axis 
position  performance  for  one  of  two  fragments  with  «1.8  ft 
spacing  between  their  trajectories,  while  the  bottom  plot  is 
the  same  except  for  ~2 . 6  ft  spacing.  The  relatively  degraded 
performance  of  the  top  plot  is  attributed  to  both  fragments' 
position  measurements  falling  within  both  filters'  validation 
gates.  This  requires  the  JPDA  algorithm  to  jointly 


This  jointly  distributed  series  of  "weak”  updates  creates 


growing  estimation  errors  along  with  filter  covariance 
divergence  that  correctly  indicates  the  filter  is  less 
confident  of  its  estimate.  In  contrast,  the  lower  plot  of 
Figure  5.5,  illustrates  the  case  where  each  fragment's 
measurement  falls  solely  within  its  own  tracking  filter's 
validation  gate  and  getting  "full  weight;"  resulting  in  low 
error  and  convergence  to  steady-state. 

Although  the  top  plot  shows  less-than-perfect  tracking, 
it  does  have  the  correct  track  for  the  correct  fragment  that 
it  started  with.  Therefore,  the  backward-running  RTS 
smoother  does  a  good  job  at  reducing  these  inaccuracies  as 
shown  in  Figure  5.6. 


Figure  5.6.  RTS  Smoother  Output  from  Top  Plot  of  Figure  5.5. 
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The  resulting  final  output  of  the  RTS  smoother  at  time  t 
may  now  be  fed  back  into  the  EKF/JPDA/RTS  package  as  initial 


conditions  and  the  whole  reduction  process  repeated.  Time 
constraints  in  this  study's  completion  prevented  the 
necessary  software  additions  to  do  such  a  "loop-around"  for 
multiple  passes;  therefore,  it  was  not  attempted. 

The  overall  review  of  the  category  two  performance 
plots.  Figures  C.55  through  C.216,  Appendix  C,  indicated  no 
significant  performance  degradation  compared  to  the  category 
one  plots,  other  than  just  discussed.  Comparing  RTS 
smoother  outputs  displayed  in  Table  5.1  with  Tables  5.2A, 
5.2B,  and  5.2C  indicates  near  identically  excellent 
performance  except  for  the  acceleration  estimates  (poor  in 
both  categories) . 


Table  V.IIA.  Category  2,  Fragment  'A',  Worse-Case  Error 
Performances 


For  initial  speed: 

lx<t0> 

=  3000. ( 

)  ft/sec 

X-Axis 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

+0.0008 

N/A 

-0.0002 

N/A 

-0.0021 

N/A 

Velocity 

+  1.25 

0.052 

-0.10 

00 

+  1.4 

0.164 

Accel . 

+  1100 

6.0 

-0.14 

00 

-800 

13.3 

For  initial  speed: 

ESS] 

=  6000. ( 

)  ft/ sec 

X-Axis 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0013 

N/A 

+0.0003 

N/A 

+0.0018 

N/A 

Velocity 

+  1.0 

0.020 

+0.25 

00 

+3.2 

0.228 

Accel . 

+20.0 

0.29 

0 

00 

+3850 

15.3 

For  initial  speed: 

|v(tc) 

=  10000 

0  ft/sec 

— 

X-Axis 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0019 

N/A 

+0.0005 

N/A 

-0.0011 

N/A 

Velocity 

+2.0 

0.024 

+0.7 

00 

+1.0 

0.036 

Accel . 

-6500 

3.22 

0 

oo 

+4000 

5.75 
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Table  V.IIB.  Category  2,  Fragment  'B',  Worse-Case  Error 
Performances 


For  initial  speed: 


X-Axis 

State  Error  Percent 


Position  +0.0005  N/A 


Velocity  -1.0  0.40 


Accel.  +600  3.3 


|v(t  ) I  =  3000.0  ft/sec 


Y-Axis 

Z-Axis 

Error 

Percent 

Error 

Percent 

+0.0014 

N/A 

+0.0006 

N/A 

+0.28 

0.636 

o 

• 

o 

i 

0.097 

+10.0 

3.25 

+  500 

9.34 

For  initial  speed:  |v(t  ) |  =  6000.0  ft/sec 


X-Axis 

Y-Axis 

Z-Axis 

Error 

Percent 

Error 

Percent 

Error 

Percent 

-0.0015 

N/A 

+0.0008 

N/A 

+0.0008 

N/A 

-3.2 

0.065 

+0.25 

0.284 

-3.0 

0.208 

-4000 

5.40 

-100 

7.74 

+4500 

21.2 

For  initial  speed:  |v(t  ) |  =  10000.0  ft/sec 


X-Axis 

Y-Axis 

Z-i 

Vxis 

Error 

Percent 

Error 

Percent 

Error 

Percent 

-0.003 

N/A 

+0.002 

N/A 

+0.0006 

N/A 

+10.0 

0.119 

+1.8 

1.26 

0.050 

-17000 

8.28 

+500 

14.6 

+6000 

10.2 

Table  V.IIC.  Category  2,  Fragment  'C',  Worse-Case  Error 
Performances 


For  initial  speed:  |v(t  )|  =  3000.0  ft/sec 


X-Axis 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

+0.0006 

N/A 

-0.0009 

N/A 

-0.0008 

N/A 

Velocity 

-0.8 

0.032 

+  0.42 

0.984 

-0.50 

0.065 

Accel . 

+500 

2.71 

o 

• 

o 

tM 

l 

6.54 

-250 

4.64 

For  initial  speed:  |v(t  ) |  =  6000.0  ft/sec 


X-Axis  Y-Axis  Z-Axis 


State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.001 

N/A 

-0.0017 

N/A 

-0.0016 

N/A 

Velocity 

0.030 

-2.5 

2.86 

+  1.7 

0.111 

Accel . 

+  1000 

1.36 

+550 

42.8 

-2900 

13.5 

For  initial  speed:  |v(tQ) |  -  20000. 0  ft /sec 


X-Axis  Y-Axis  Z-Axis 


State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0018 

N/A 

-0.0008 

N/A 

-0.0014 

N/A 

Velocity 

-5.0 

0.061 

+0.5 

0.347 

+2.1 

0.083 

Accel . 

-6000 

2.94 

-100 

2.92 

-5000 

8.40 

One  other  impact  of  tracking  multiple  fragments  versus  a 
single  fragment  was  the  increase  in  computer  run  time.  For  a 
single  fragment  run,  the  EKF/JPDA/RTS  package  needs  only  to 
run  one  EKF  and  the  measurement  validation  process  simplifies 
to  one  g-sigma  test.  Such  simulation  runs  were  completed  in 
a  matter  of  minutes  on  a  VAX  8650  host.  In  contrast,  for 
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tracking  three  fragments,  the  package  must  run  nine 
simultaneous  EKFs  and  test  each  filter's  expected 
measurement,  h[*,*],  against  nine  candidate  measurements,  z^. 
Following  the  event  algebra  discussion  in  Appendix  B,  for  a 
no-clutter  scenario,  Eq. (B.lb)  defines  the  JPDA  algorithm 
must  "step  through"  over  9!  =  362,880  event  permutations. 
However,  for  this  study,  which  includes  clutter 
consideration,  Eq. (B.7) ,  Appendix  B,  increases  the  total 
events  to  17,572,114  permutation/combination  events  for 
each  update  cycle.  The  resulting  computational  time  commonly 
extended  to  over  24  hours.  For  a  case  where  four  valid 
fragments  are  to  be  tracked,  Eq. (B.7)  defines  a  total  of 
6,199,668,952,530,000  possible  events  of  measurement-to- 
target  permutations  and  all  clutter  combinations.  Applying 
this  to  a  hypothetical  computer  that  could  do  ten-million 
permutations  per  second,  the  above  run  would  require  over  19 
years  to  do  one  update.  Needless  to  say,  the  event 
generation  routine  developed  in  this  study  that  must  "step 
through"  every  possible  event  to  find  all  Q(X)  event 
matrices  for  a  given  Q  validation  matrix  must  be  optimized 
by  some  manner  for  practical  data  reduction  use. 

V.2.C  Category  Three .  This  category,  with  its  curving 
and  crossing  trajectories,  combined  with  the  presence  of  non- 
trackable  clutter  fragments,  posed  the  "worse  case"  challenge 
to  the  EKF/JPDA/RTS  package.  It  is  emphasized  here,  that 
these  category  three  runs  were  modelled  with  unrealistically 
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large  magnitude  lateral  accelerations  in  order  to  get  the 
trajectories  to  cross.  This  was  intended  to  thoroughly  test 
the  EKF/JPDA/RTS  package's  tracking  capabilities  under 
adverse  conditions.  Such  curving  acceleration  magnitudes  are 
highly  unlikely  in  a  true  test  blast  fragment  dispersion. 

The  resulting  runs  (as  expected)  gave  the  least  overall 
performance  of  the  three  categories.  For  these  runs,  two 
fragments  were  given  a  substantially  large  lateral 
acceleration  (mostly  in  the  Y-axis) ,  such  that  their 
trajectories  curved  and  crossed  tracks  approximately  one  foot 
from  each  other,  while  within  both  cameras'  fields-of-view. 
Simultaneously,  two  clutter  fragments,  each  only  visible  to 
one  camera,  were  generated  to  create  "extra"  measurements 
that  the  JPDA  algorithm  had  to  process  out.  Specifically, 
the  clutter  fragment  trajectories  were  designed  such  that  one 
passed  through  the  far-field  view  of  camera  #1  and  the  other 
passed  through  the  near-field  view  of  camera  #2,  as  shown  in 
Figure  5.7.  Such  a  clutter  situation  is  quite  realistic  for 
the  arena  test  application. 

From  Figures  C.217  through  C.324  in  Appendix  C,  one  can 
quickly  determine  the  EKFs  experiencing  severe  transient 
conditions  since  they  where  initialized  assuming  straight- 
line  trajectories  from  the  arena  center. 


Figure  5.7.  Clutter  Fragment  Trajectory  Locations. 


For  the  slow,  3000  ft/sec  fragment  run,  the  EKFs  receive 
enough  measurements  to  recover  good  tracking  in  the  Xj  and 
Zj  axis  directions  since  their  lateral  acceleration  and 
velocity  components  are  minimal.  However,  the  defined  arena 


geometry  combined  with  the  "straight  path"  assumption  leads 
the  EKF  initialization  procedure  discussed  in  Chapter  3  to 


expect  the  greatest  fragment  acceleration  in  the  Xj  axis 
(radial)  direction  and  near-zero  acceleration  in  the  Yj 
axis  (tangential)  direction  (review  Figures  A. 3,  A. 4  and  A. 5 
in  Appendix  A) .  Therefore,  the  unexpected  large  Yj 
acceleration  "fools"  the  EKFs  as  they  propagate  a  straight 
trajectory  while  the  fragment  is  actually  curving. 


Fortunately,  the  a^^  threshold  feature  of  the 
initialization  process  prevents  the  Yj  acceleration 
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component  from  being  set  lower  than  one-third  the  Xj 
acceleration  value.  This  keeps  the  Kalman  gains  for  the  Yj 
axis  direction  high  enough  to  help  "pull"  the  EKF  states 
toward  the  correct  values  from  the  succeeding  measurement 
updates.  Although  the  presented  tracks  in  the  Yj  axis 
direction  contain  diverging  errors  that  surpass  the  EKFs ' 
standard  deviation  values,  the  results  obtained  from  the  RTS 
smoother  outputs  greatly  reduced  these  inaccuracies.  These 
RTS  outputs  for  time  t  may  now  be  substituted  as  initial 
conditions  back  into  the  EKFs  for  reprocessing,  where  this 
time,  the  EKFs  are  better  initialized  to  reduce  the  previous 
transient  behavior.  Such  a  multi-cycle  process  seems 
especially  applicable  to  the  6000  and  10000  ft/sec  fragment 
runs  since  they  need  the  "extra"  recursion  cycles  to 
approach  steady-state. 

Table  5.3A  and  5.3B  summarize  the  category  three  results 
where  the  Yj  axis  shows  the  only  significant  degradation  in 
performance.  Similar  to  categories  one  and  two,  the 
acceleration  estimates  for  this  category  are  considered 
unreliable,  and  are  extremely  poor  in  the  Yj  axis 
direction. 

Special  notice  is  brought  to  the  fact  that  no  degrading 
effects  could  be  found  from  the  presence  of  the  two  clutter 
fragments  (except  computational  run-time) .  This  makes 
perfect  sense  since  candidate  measurements  created  by  the 
clutter  should  consistently  "fail"  all  validation  tests. 


For  initial  speed:  |v(t  ) |  =  3000.0  ft/sec 


X-Axis 


Y-Axis 


Z-Axis 


State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.002 

N/A 

+0.009 

N/A 

+0.0003 

N/A 

Velocity 

+4.5 

0.18 

+  11.0 

mm am 

0.91 

Accel . 

-2150 

10.3 

+  5500 

99.5 

-700 

12.5 

For  initial  speed: 

lY(tc) 

=  6000.0  ft/sec 

X-i 

^xis 

y -; 

Ixis 

2-1 

^xis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0035 

N/A 

+0.0102 

N/A 

+0.0007 

N/A 

Velocity 

+8.2 

0.16 

+27.0 

21.3 

-2.8 

0.18 

Accel . 

+9300 

14.4 

+22000 

99.5 

+3300 

14.1 

For  initial  speed: 

ls<V 

=  10000.0  ft/sec 

X-Axis 

Y-J 

^xis 

2-1 

^xis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0058 

N/A 

-0.025 

N/A 

-0.0024 

N/A 

Velocity 

-33.0 

0.04 

-70.0 

34.6 

+  1.9 

0.08 

Accel . 

+56000 

24.2 

-61000 

99.7 

-5000 

9.20 

Table  V.IIIB.  Category  3,  Fragment  ' B ' ,  Worse-Case  Error 
Performances 

For  initial  speed:  |v(tQ) |  =  3000.0  ft/sec 

X-Axis  Y-Axis  Z-Axis 


State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

-0.0006 

N/A 

-0.007 

N/A 

-0.0005 

N/A 

Velocity 

-1.0 

0.042 

-10.5 

16.4 

o 

01 

• 

o 

1 

0.12 

Accel . 

-700 

3.50 

-5000 

90.8 

+720 

13.9 

For  initial  speed:  jv(t  ) |  =  6000.0  ft/sec 


X-Axis  Y-Axis  Z-Axis 


State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

+0.0045 

N/A 

-0.0085 

N/A 

-0.0011 

N/A 

Velocity 

-9.0 

0.18 

-26.0 

20.7 

-1.8 

0.13 

Accel . 

+8800 

10.5 

-22000 

99.8 

-2000 

10.2 

For  initial  speed:  jv(t  ) |  =  10000.0  ft/sec 


X-Axis 

Y-Axis 

Z-Axis 

State 

Error 

Percent 

Error 

Percent 

Error 

Percent 

Position 

+0.008 

N/A 

+0.027 

N/A 

+0.003 

N/A 

Velocity 

+41.0 

0.48 

+72.0 

35.3 

+4.7 

0.19 

Accel . 

-61000 

34.3 

+62000 

100 

+8000 

12.4 

Therefore,  when  a  JPDA  event  matrix  matches  a  clutter 
measurement,  the  algorithm  will  simply  step  past  it  and  go  on 
to  the  next  event.  This  precludes  a  "weight"  being  assigned 
to  the  clutter  measurement  and  subsequent  "weakening"  of  the 
resulting  update. 
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This  chapter  has  presented  the  three  categories  of 
fragment  trajectories  simulated  in  this  study  along  with  a 
special  run  to  indicate  the  effects  of  20  versus  100  Monte- 
Carlo  runs  on  simulation  results.  For  each  category, 
specific  performance  highlights  and  shortcomings  were 
discussed  along  with  any  "attempted  fixes."  Further 
discussion  described  the  computational  workload  shortfall 
encountered  with  this  study.  One  other  special  run  included 
in  this  chapter  compared  the  tracking  performance  of  "close" 
versus  "distant"  fragment  trajectories,  relative  to  each 
other . 


This  study  has  developed,  and  simulated  a  multiple 
fragment  tracking  system  for  use  in  munitions  testing. 


Simulation  analysis  has  shown  overall  positive  results  and 
the  usefulness  of  the  JPDA  algorithm.  Some  important 
observations  learned  from  the  research,  design,  simulation, 
and  analysis  are  given  below: 

VI . 1  Conclusions 

1. )  The  designated  design  specification  for  a  CCD 
camera  of  512  x  512  pixel  resolution,  operating  at  5000 
frames  per  second,  is  a  very  demanding  specification  to  meet 
because  of  the  multi-gigahertz  (GHz)  signal  processing 
speeds  required  to  read,  convert  (analog  to  digital) ,  and 
store  the  image  data  from  the  camera  chip  to  a  suitable 
storage  medium.  Present  technology,  non-imaging  CCD  devices 
have  been  operated  at  the  GHz  range,  but  no  report 
literature  was  found  to  support  the  5000  frame/sec  rate  as 
truly  achievable  with  present  technologies. 

2. )  The  orthogonal  sensor  (camera)  geometry  simplified 
the  applied  measurement  model  by  eliminating  cross-covariance 
noise  terms. 

3. )  The  linear,  first-order,  Gauss-Markov  acceleration 
model  employed  within  the  EKFs  proved  to  be  a  very  close 
match  to  the  nonlinear  fragment  acceleration  truth  model, 

Eq. (3.1) . 
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4. )  Tracking  error  performance  of  straight-trajectory 
fragments  in  position  often  out-performed  the  camera  pixel 
resolution  by  an  order  of  magnitude,  even  though  the  camera 
measurements  were  corrupted  by  a  three-sigma  noise  level  of 
one  pixel.  Similarly,  velocity  tracking  errors  rarely 
surpassed  0.2%  of  the  true  state  value. 

5. )  The  acceleration  state  estimates  for  all  simulation 
runs  proved  inconsistent  and  unreliable  to  the  point  that 
they  were  determined  unusable. 

6. )  The  EKF  initialization  process  developed  in  this 
study  performed  extremely  will  for  straight-path  and  mildly 
curving  fragments.  This  process  also  demonstrated  good 
matching  of  the  acceleration  time  constant  parameters,  r  , 

t  ,  and  r z,  to  the  true  acceleration  decay  rates. 

7. )  The  slow  adaptation  rate  of  the  EKFs  resulted  in 
their  inability  to  hold  track  on  strong-curving  fragments. 

8. )  The  inclusion  of  the  RTS  smoother  greatly  improved 
the  quality  of  the  estimated  outputs  without  significantly 
increasing  the  computational  run-time. 

9. )  Although  this  study  developed  an  operational 
software  algorithm  to  generate  all  possible  events  for  a 
given  validation  matrix,  its  computational  load  made  it 
impractical  for  applications  involving  more  than  three 
fragments  visible  in  each  camera  (nine  measurement 
combinations) . 
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10. )  The  JPDA  validation  gate  process  employing  the  "g- 
sigma"  test  proved  very  effective  at  eliminating  undesirable 
candidate  measurements,  such  as  the  simulated  clutter 
fragments . 

11. )  The  degraded  performance  of  the  JPDA  algorithm  for 
"close"  trajectories  could  not  be  improved  by  simply 
adjusting  the  validation  gate  size,  g.  The  gate  size  could 
be  reduced  to  exclude  the  "other  track"  measurements,  but  at 
the  cost  of  the  gate  being  so  small,  that  the  measurement 
noise  would  often  "jig"  the  correct  track  measurements 
outside  the  proper  gate.  This  resulted  in  a  significantly 
reduced  number  of  updates  and  complete  loss  of  track. 

VI. 2  Recommendations 

Many  of  the  assumptions  made  while  completing  this  study 
were  made  either  to  simplify  the  problem  or  because  no  data 
could  be  found  to  support  alternative  choices.  The  following 
recommendations  consider  variations  to  the  general 
assumptions  outline  in  Chapter  I: 

1. )  The  simplified,  spherical  shaped,  fragment  drag 
model  needs  to  be  expanded  to  non-spherical  shapes,  such  as 
platelets,  wedges,  etc. 

2. )  The  simplified,  white,  Gaussian  measurement  noise 
needs  to  expanded  to  include  biases  for  camera  misalignment 
and  time-correlated  noise  components  for  camera  vibration, 
etc. 
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3. )  The  image  pre-processing,  that  was  assumed  already 
present,  must  be  designed  that  will  extract  the  pertinent 
fragment  data  from  the  raw  digitized  image  data.  This 
process  must  extract  as  a  minimum,  the  quantity  of  fragments 
visible,  fragment  diameters,  and  fragment  image  centroid 
coordinates. 

4. )  This  study  set  the  JPDA  clutter  density  value,  C, 
to  unity  for  all  simulation  runs  since  no  data  was  available 
to  do  otherwise.  Likewise,  the  camera  probabilities  of 
detecting  fragments  within  their  f ield-of-view,  PD,  were 
arbitrarily  set  at  0.95  simply  because  it  was  a  "nice” 
number.  Variations  to  these  parameters  need  investigation 
and  their  resulting  effects  on  the  JPDA  algorithm 
performance. 

5. )  This  study,  through  ad  hoc  tuning,  determined  that 
validation  gate  values,  g,  were  best  set  at  3.0  for  straight- 
trajectory  fragments  and  4.0  for  curving  trajectory 
fragments.  A  more  analytic  means  of  determining  this 
parameter  as  a  function  measurement  noise  strength,  R(*),  and 
other  input  sources  should  be  investigated. 

6. )  This  study  held  the  acceleration  time  constants, 

t  ,  t  ,  and  r  ,  constant  once  they  were  initialized  at  t  . 
x  y  z  o 

The  consideration  to  augment  these  to  the  EKF  as  estimated 
parameters  may  be  considered. 

7. )  The  simulation  software  developed  in  this  study 
should  be  modified  to  allow  the  RTS  smoother  output  at  time 


A 


\ 

n 
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t  be  "looped  back"  into  the  EKF  input  as  initial  conditions. 
This  would  facilitate  multiple,  bi-directional  processing 
runs  for  a  given  set  of  measurements. 

8.)  The  algorithm  developed  in  this  study  to  generate 
each  possible  event  matrix,  n(X) ,  for  a  given  validation 
matrix,  n,  must  be  rewritten  to  be  less  computationally 
burdensome . 
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APPENDIX  A:  Measurement  Geometry 

This  appendix  details  the  development  of  the  applied 
observation  and  measurement  geometry,  including  specifics 
defining  reference  frames  and  transform  functions  between 
reference  frames. 

The  general  arena  setup  involves  the  following  set  of 
three  reference  frames: 

Frame  1)  A  two  dimensional  (2D)  x-y  coordinate  plane 
for  each  camera  where  the  origin  is  located  at  the  center  of 
the  image  plane  (camera  bore-sight) : 


Figure  Al.  Camera  Object  Plane  Coordinate  Frame 


For  the  above  image  coordinates,  the  x  and  y  designators 
are  subscripted  by  the  camera  number.  For  example  camera  one 
would  have  a  coordinate  pair  listed  as  (x1#  y1)  and  camera 
two  would  have  its  coordinate  pair  listed  as  (x2,  y2) • 

The  width  and  height  of  the  image  plane  (modelled  to  be 
square,  so  width  =  height)  is  determined  by  the  focal 
distance,  D,  and  the  camera's  angular  FOV,  a,  determined  by 
camera  lens  selection  as  depicted  in  Figure  A2: 
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Figure  A2.  Camera  FOV  Geometry 

The  function  relating  object  plane  Size  to  D  and  a  in 
Figure  A2: 

Size  =  2  D  TAN(a/2)  (A-l) 

From  this  equation,  various  image  plane  sizes  for 
various  resolutions  or  arena  set-up  requirements  may  be  met 
by  changing  camera  lenses  or  adjusting  D  in  the  arena  set¬ 
up. 

Frame  2)  Each  camera  pair  (x^  y^,  x2,  y2)  defines  a 
three  dimensional  (3D)  right-hand  intermediate  reference 
frame  (Xj,  Yj,  Zj)  such  that  the  object  planes  of  the  two 
cameras  orthogonally  intersect  and  both  camera  object  frame 
origins  and  the  intermediate  frame  origin  are  at  the  same 
point.  The  Zj  direction  is  pointed  in  the  local  vertical 
(up) .  The  Xj  direction  is  pointed  radially  outward  from 
the  arena  center.  The  remaining  Yj  direction  is  set  by 
the  right-hand  rule,  Zj  cross  Xj  equals  Yj.  See  Figures 
A3 ,  A4  and  A5 . 
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Figure  A3.  Overhead  View  of  Arena  Geometry 


Figure  A4.  Radial  View  (looking  inward)  of  Arena  Geometry 


The  transformation  from  the  2D  camera  image  coordinates 
(C-frame)  to  the  3D  intermediate  frame  (I-frame)  requires  a 
set  of  three  nonlinear  functions.  The  derivation  of  these 
functions  is  based  on  combining  the  Pythagorean  theorem, 
properties  of  similar  triangles,  intersection  point  of  two 
lines,  and  numerous  algebraic  substitutions  to  yield: 

2  D  x.  x_ 

X  =  - -  (A-2) 

'  xi  y2  ‘  x2  yi  +  D  x2  _  D  xi 


Y 


I 


p2  (  yi  -  y2  ) 

R  (  y1  y2  -  D2  ) 


(A-3) 


J 


fl 


I 


i 


fe 


i 


i 
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(A-4 ) 


ZI  = 


where 


D  [  y2  (2y1+P)+Py1] 

R  (  D2  -  yx  y2  ) 

1/2 


R  =  (  2.0  ) 


Note  that  Eq. (A-2)  is  undefined  when  x1 ,  y  ,  x2,  and  y2 
are  all  zero.  This  condition  must  be  checked  for  in  any 
software  implementations. 

To  transform  from  the  I-frame  back  to  the  C-frame 
requires  the  above  system  of  equations  be  algebraically 
manipulated  and  solved  for  x^,  y1#  x2  and  y2  (Note:  There  is 
a  second  solution  for  Yj  other  than  that  shown  as  Eq. (A-3) ; 
therefore,  four  equations  do  exist  to  solve  for  the  four 
unknowns) .  The  resulting  inverse  transform  equations  are: 


X1  " 


2  D  Xj  (  R  Yj  +  D  ) 


2  (  Yj  Zj  +  Yj  +  D2  )  +  R  D  (  3  Yj  +  Zx  ) 


*1  = 


R  D  (  Zj  -  Yj  ) 

R  (  Zj  +  Yj  )  +  2  D 


y2 


R  D  Xj  (  zx  +  Yt  ) 

R  (  Zj  ■  Yj  )  +  2  D 


D  X 


I 


D  -  R  Yj 

R  D  (  ZT  4  YT  ) 

R  (  Zj  -  Yj  )  +  2  D 


(A-5) 


(A-6) 


(A-7) 


(A-8) 
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Frame  3)  A  right-hand  3D  main  reference  frame  or  world 
reference  frame  Yw,  Zw  with  its  origin  at  the  center  of 

the  arena  (also  location  of  test  munition)  (again  review 
Figures  A. 3,  A. 4,  and  A. 5).  The  Zw  axis  is  pointed  to 
local  vertical  (up,  parallel  to  Zj)  .  The  XW~YW  plane  is 
therefore  horizontal.  The  direction  to  which  points  is 

arbitrary  (such  as,  pointing  north)  since  this  W-frame  is  the 
master  "inertial”  frame  to  which  all  previously  defined 
reference  frames  eventually  transform. 

The  transformation  function  from  the  I-frame  to  the  W- 
frame  simply  involves  a  single  direction  cosine  matrix  (DCM) 
operation  for  fragment  velocity  and  acceleration  vectors, 
while  position  vectors  also  require  an  additional 
translation  vector  component.  As  indicated  in  the  previous 
reference  frame  definitions,  the  Zw  and  Z1  directions  are 
parallel  making  them  the  "pivot"  for  the  rotation  process. 

Based  on  the  selected  location  of  each  camera  pair,  the 
correspondingly  created  I-frame  origin  is  located  at  angle  0 
from  the  X^  axis  following  standard  right-hand  rotation 


(see  Figure  A3) .  The  desired  DCM,  H 


,W 


becomes : 


cos  e  -sin  e 
sin  e  cos  e 


o 

o 

l 


(A-9) 


and  is  applied  as: 
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(A-10) 


YW  =  Mj  Y1 


where  the  inverse  transform  from  the  W-frame  to  the  I-frame 
is  the  transpose  of  the  above  DCM  to  give  : 


cos  e  sin  e  o 

-sin  e  cos  e  o 

001 


(A-ll) 


such  that: 


**  -  K  *w 


(A-12) 


When  transforming  fragment  position  vectors,  a  position 
translation  vector  from  the  W-frame  origin  to  the  I-frame 
origin  must  also  be  included.  The  components  of  this 
translation  vector,  t  are  the  x,  y  and  z  projections  of 
this  translation  vector  onto  the  Xw,  Yw,  and  Zw  axes 
respectively: 


Let:  Rd  =  Arena  radius  (from  arena  center  to  camera  lens 

measured  along  ground  plane)  (Figure  A. 3). 

Sz  =  The  vertical  displacement  of  the  W-frame  and 

I-frame  origins  (preselected  for  particular  arena 
set-up)  (Figure  A. 5). 


(  Rd2  -  Sz2 

j1/2  cos  e 

(A-13) 

(  Rd2  -  $z2 

)1/2  sin  e 

(A-14 ) 

Sz 

(A-15) 
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When  transforming  an  I-frame  position  vector,  p1,  to  the 
w 

W-frame,  p  ,  the  t  vector  is  summed  to  the  DCM* I-frame 
product  as: 

W  W  T 

p=  M"  p  +  t  (A-16) 

When  transforming  from  W-frame  to  I-frame  the  t  vector  is 
subtracted  from  the  DCM* W-frame  product  as: 

p1  =  pW  -  t  (A-17) 

This  concludes  the  reference  frames  and  associated 
transforms  discussion  pertinent  to  this  research  effort. 
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APPENDIX  B:  Event  Algebra 

The  following  derivations  develop  the  worst-case  I 

situation  for  a  given  validation  matrix,  Q,  that  is  fully 
unit  populated  (all  ones) .  Although  such  an  occurrence  is 
highly  unlikely  the  JPDA  algorithm  must  "assume  the  worse"  I 

since  there  is  no  "a  priori"  information  regarding  which  or 
how  many  of  the  candidate  measurement-to-target  associations 
are  true.  I 


Begin  with  the  following  variable  definitions: 

T  =  The  number  of  candidate  non-clutter  targets  (tracks) .  \ 

M  =  The  number  of  candidate  measurements, 
f  =  The  number  of  false  measurements  (clutter) . 

E  =  The  integer  number  of  possible  feasible  events.  I 

t  =  The  index  counter  for  candidate  tracks;  0  <  t  <  T  . 
j  =  The  index  counter  for  candidate  measurements;  1  <  j  <  M  . 


where  also  the  number  of  permutations,  P,  of  n  objects 
taken  r  at  a  time  may  be  defined  as: 


n! 


nPr 


(n  -  r) ! 


(B.l) 


and  the  number  of  combinations,  C,  of  n  objects  taken  r 
at  a  time  may  be  defined  as: 


n! 


C  = 
n  r 


rl  (n  -  r)  1 


(B.2) 


Recall  from  Chapter  2,  a  feasible  event  is  defined  where  no 
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more  than  one  measurement  may  originate  from  each  target 
(24:176) : 


j  ^  r  and  t^  >  0  implies  tj  ^  tr  (B.3 

Also  repeating  the  rules  defining  an  event,  Q(X) ,  from  the 
given  validation  matrix,  Q  : 


1)  Scan  Q  by  rows  and  pick  one  unit  per  row  for 
Q(X). 

2)  Only  one  unit  from  each  column  i  >  1  can  be  taken 
(i.e.  at  most  one  measurement  could  have  originated 
from  a  target) .  The  number  of  units  from  column 

i  =  0  is  not  restricted. 


From  these  definitions  and  beginning  with  a  simple  "no 
clutter"  case,  the  number  of  targets,  T,  and  the  number  of 
measurements,  M,  replace  n  and  r  respectively  in 
Eq. (B.l)  to  give  an  fi  matrix  and  the  total  number  of 
events,  E,  as: 


T! 

(T  -  M)  ! 
for  M  <  T 


(B.la 


where  two  classifications  describe  as: 

Case  1:  Eq.(B.la)  applies  directly  for  no-clutter  and 
M  <  T  . 

Case  2:  Same  as  Case  1  except  when  M  =  T  ,  Eq.(B.la) 
simplifies  to: 


T 


tEt  =  T!  (B. lb) 

For  the  situation  where  M  >  T,  the  feasibility  criteria 
is  violated  unless  provisions  are  made  to  account  for 
clutter.  This  is  done  by  augmenting  an  extra  clutter,  tQ 
column  to  the  front  of  the  validation  matrix  as: 


Eq.(B.la)  is  modified  to  include  the  additional  clutter 
events  where  M  is  replaced  with  (M  -  f)  to  include 
assumed  feasible  measurements  minus  clutter  measurements: 


T! 

tEm  f  =  -  ;  for  M  >  T,  M  <  T  or  M  =  T  (B.4) 

'  (T  -  M  +  f )  ! 

where  for  cases  of  M  >  T  then  f  must  be  within  the  range: 

(M  -  T)  <  f  <  M 

However,  Eq.(B.4)  accounts  only  for  the  feasible 
measurements  while  neglecting  the  clutter  possibilities.  The 
additional  combinations  of  clutter  events  are  accounted  for 
by  inserting  M  and  f  in  place  of  n  and  r 


4 
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respectively  in  Eq.(B.2)  to  give: 


Mi 


c  = 

vrf 


f  i  (M  -  f )  ! 


(B.  5) 


The  product  of  Eqs.(B.4)  and  (B.5)  results  in  the  total 
number  of  possible  events  for  any  ratio  of  T  to  M  and  a 
fixed  number,  f,  of  clutter  measurements: 


TEM,f 


Mi  Ti 

f !  (M  -  f)  !  (T  -  M  +  f  )  1 


(B.  6) 


where  f  must  be  in  the  range: 


(M  -  T)  <  f  <  M 

Note  that  Eq. (B. 6)  reduces  to  Eq.(B.la)  when  f  =  0  (no 
clutter)  and  Eq.(B.lb)  when  f  =  0  and  T  =  M  . 

Since  the  number  of  clutter  measurements  is  not  fixed, 
all  possible  numbers  of  clutter  measurements  must  be 
accounted  for  where  the  number  of  possible  clutter 
measurements  may  fall  in  the  range:  (M  -  T)  <  f  <  M  . 
Therefore,  the  absolute  total  number  of  possible  events 
including  clutter  for  any  given  T  to  M  ratio  results  as  a 
summation  of  multiple  cases  of  Eq.(B.6)  over  the  range  of 
clutter  possibilities: 


Total 


M 

Mi  T!  V 
q  =  0 


1 


f !  (M  -  f )  !  (T  -  M  +  f )  ! 
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where 


q  =  0  for  (M  -  T)  <  0 

q  =  (M  -  T)  for  (M  -  T)  >  0 


(B.  7 ) 


The  final  quantity,  E,  given  in  Eq.(B.7)  grows  very 
rapidly  for  increasing  numbers  of  targets  and  measurements 
due  to  the  factorial  operators.  Re-emphasizing  that  a  fully 
populated  validation  matrix  is  extremely  unlikely,  the  JPDA 
event  generator  algorithm  designed  in  this  study,  must  still 
"step  through"  each  possible  permutation-combination  event, 
Q(X) ,  and  test  it  for  a  match  against  the  present  validation 
matrix,  n.  The  resulting  computational  workload  may  become 
overwhelming. 


Simulation  Output  Plots 


The  plots  contained  m  this  appendix  were  generated 
using  the  MATRIX^  interactive  software  package  plotting 
utility  (32).  The  following  lists  summarize  the  simulation 
parameters  that  all  the  contained  plots  share  in  common: 


1. )  Created  from  20  Monte-Carlo  runs. 

2. )  Object  plane  size  and  resolution  level: 

10.0  x  10.0  ft  @  512  x  512  pixels 
(»0.02  x  «0. 02  ft  per  pixel  side) 

3 .  )  EKF  measurement  model  expected  three-sigma  noise 

level : 

0.02  ft 

4. )  JPDA  sensor  probability  of  detection:  PQ  =  0.95 

5. )  JPDA  Poisson  clutter  density:  c  *  1.0  ft-4 

6. )  True  position  and  diameter  measurement  three-sigma 

noise  level: 

0.02  ft 


The  nomenclature  used  the  label  the  vertical  axis  of 
each  plot  is  derived  from  the  following  list: 


XPERR 


XVERR 


XAERR 


YPERR 


YVERR 


YAERR 


The  X-axis  position  error  performance. 

The  X-axis  velocity  error  performance. 

The  X-axis  acceleration  error  performance 
The  Y-axis  position  error  performance. 

The  Y-axis  velocity  error  performance. 

The  Y-axis  acceleration  error  performance 


ZPERR  :  The  Z-axis  position  error  performance. 

ZVERR  :  The  Z-axis  velocity  error  performance. 

ZAERR  :  The  Z-axis  acceleration  error  performance. 

For  all  plots,  the  horizontal  axis  labeled  SEC  is  the 
time  in  seconds  following  detonation.  Other  parameters,  such 
as  fragment  initial  speed,  trajectory  category,  and  fragment 
index  label  A,  B,  or  C  are  indicated  on  each  plot. 
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Figure  C.  3.  Tracking  Error  Plot,  Category  1,  v(to)  -  3000  ft/sec. 


mdk 


Figure  C.  8.  Tracking  Error  Plot,  Category  1,  v(to)  -  3000  ft/sec. 
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C.  10.  Tracking  Error  Plot,  category  1,  v(to) 
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Figure  C 


13.  Tracking  Error  Plot, 
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Figure  C.  16. 


Tracking  Error  Plot,  Category  1,  v(to) 


3000  ft/sec. 
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Figure  c.  26.  Tracking  Error  Plot,  Category  1,  v(to)  *  6000  ft/sec. 
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Figure  C.  35 


Tracking  Error  Plot,  Category  1,  v(to) 
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C.  42.  Tracking  Error  Plot,  Category  1,  v(to)  ■  10000  ft/sec 
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Figure  C.  43.  Tracking  Error  Plot,  Category  1,  v(to)  *  10000  ft/sec 
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Figure  C.  46.  Tracking  Error  Plo 
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Figure  c.  «a.  Tracking  Error  Plot,  Category  1,  v(to)  -  loooo  ft/eec 
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Figure  C.  49.  Tracking  Error  Plot,  Category  1,  v(to) 
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Figure  C.  55.  Tracking  Error  Plot,  Category  2,  v(to) 


3000  ft/sec 


Figure  C.  57.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.  59.  Tracking  Error  Plot,  Category  2,  v(to)  «  3000  ft/sec 
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Figure  C.  60.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.  61.  Tracking  Error  Plot,  Category  2,  v(to) 


3000  ft/sec 
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Figure  C.  67.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.  68. 


Tracking  Error  Plot,  Category  2,  v(to) 


3000  ft/sec. 
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.  85.  Tracking  Error  Plot,  Category  2,  v(to)  -  3000  ft/sec. 
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Figure  C.  90.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.108.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.110. 


Tracking  Error  Plot,  Category  2,  v(to) 


6000  ft/sec. 
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Figure  C.112.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.114.  Tracking  Error  Plot,  Category  2,  v(to) 
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C. 115.  Tracking  Error  Plot,  Category  2,  v(to)  -  6000  ft/sec. 


233 


4000 


238 


243 


dd3AX 


Figure  C.130.  Tracking  Error  Plot,  Category  2,  v(to) 


6000  ft/sec. 
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Figure  C.132.  Tracking  Error  Plot,  Category  2,  v(to)  «  6000  ft/sec 


I 


M3dA 


Figure  C.134.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.142. 


Tracking  Error  Plot,  Category  2,  v(to) 


6000  ft/sec. 
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Figure  C.144.  Tracking  Error  Plot,  Category  2,  v(to) 
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.147.  Tracking  Error  Plot,  Category  2,  v(to)  -  6000  ft/sec. 
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Figure  C.15P. 


Tracking  Error  Plot, 


Category  2,  v(to) 


6000  ft/sec. 
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Figure  C.166.  Tracking  Error  Plot,  Category  2,  v(to)  -  10000  ft/sec 
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Figure  C.180.  Tracking  Error  Plot,  Category  2,  v(to) 
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Figure  C.181.  Tracking  Error  Plot,  Category  2,  v(to)  -  10000  ft/ sec. 
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Figure  C.182.  Tracking  Error  Plot,  Category  2,  v(to)  «  10000  ft/sec. 
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Figure  C.184 


Tracking  Error  Plot,  Category  2,  v(to)  -  10000  ft/sec 
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Figure  C.186.  Tracking  Error  Plot,  Category  2,  v(to)  -  10000  ft/sec 
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C.187.  Tracking  Error  Plot,  Category  2,  v(to)  »  10000  ft/sec. 
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Figure  C.192. 


Tracking  Error  Plot,  Category  2,  v(to)  «  10000  ft/sec 
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Figure  C.195.  Tracking  Error  Plot,  Category  2,  v(to)  «  10000  ft/sec 
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Figure  C.198.  Tracking  Error  Plot,  Category  2,  v(to)  «  10000  ft/sec 


319 


I 

1 

!~ 

32 


0009 


Figure  C.224.  Tracking  Error  Plot,  Category  3,  v(to) 
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Figure  C.228.  Tracking  Error  Plot,  Category  3,  v(to) 
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Figure  C.234.  Tracking  Error  Plot,  Category  3,  v(to) 
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Figure  C.235.  Tracking  Error  Plot,  Category  3,  v(to) 
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Figure  C.239.  Tracking  Error  Plot,  Category  3,  v(to) 
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Figure  C.260.  Tracking  Error  Plot,  Category  3,  v(to)  =  6000  ft/sec. 


378 


Figure  C.264.  Tracking  Error  Plot,  Category  3,  v(to) 
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Figure  C.282.  Tracking  Error  Plot,  Category  3,  v(to)  =  6000  ft/sec 
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C. 283 .  Tracking  Error  Plot,  Category  3,  v(to)  =  6000  ft/sec. 
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Figure  C.288.  Tracking  Error  Plot,  Category  3,  v(to) 
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The  purpose  of  this  study  was  to  research  and  apply 
current  technology  electronic  data  acquisition  and  tracking 
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Previously  applied  technologies  for  munitions  testing  have 
proven  themselves  expensive  in  materials,  labor,  and  time. 

The  specific  scope  of  this  study  was  to  research  methods  to 
electronically  acquire  and  track  the  position,  velocity,  and 
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from  the  test-blast  center.  A  design  specification  for  a 
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The  theoretical  application  of  xenon  strobe  illuminated  (2.0 
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dimensional  image  measurements  at  a  2.0  microsecond  exposure, 
5000  frame/sec  rate.  The  acquired  and  assumed  noisy  fragment 
position  measurements  (recorded  digitally)  are  post-mission 
processed  through  an  Extended  Kalman  Filter  (EKF)  based  Joint 
Probability  Data  Association  (JPDA)  multiple  target/tracker 
state  estimator  followed  by  a  backward  time  Rauch-Tung- 
Striebel  (RTS)  fixed-interval  optimal  smoother.  Strong 
emphasis  was  placed  on  Monte-Carlo  computer  simulation 
testing  of  this  EKF/JPDA/RTS  tracker-smoother  algorithm. 
Representative  trajectories  of  straight,  curving,  and 
crossing  spherical  fragments  at  3000,  6000  and  10000  ft/sec 
were  modeled  and  tracked  with  promising  accuracies  in 
position  and  velocity^  The  presented  fragment  data 
acquisition  system  was  deemed  realizable  and  practical  with 
existing  technologies,  although  the  CCD  camera  5000  frame/sec 
requirement  was  found  difficult  to  obtain  reliably.  The 
initial  proposed  system  hardware  cost  will  be  high;  however, 
critical  system  components  (such  as  cameras)  survive  the  test 
blast  and  are  continuously  reusable  to  keep  overall  long-term 
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